169 USE iso_fortran_env
, only: int8
196 real:: safe_ratio=2.1
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
211 integer:: check_interval=10
212 real:: min_dx=0.0, max_dx=1e32
213 real:: refine_from_t=0.0
214 real:: refine_until_t=1e30
217 real:: max_shock=0.0, max_contact=0.0
218 real:: min_shock=0.0, min_contact=0.0
219 real:: min_compress=0.0, max_compress=0.0
220 real:: min_vorticity=0.0, max_vorticity=0.0
221 integer,
parameter:: ngradvar=10
222 real,
dimension(3):: region_llc=0.0
223 real,
dimension(3):: region_urc=0.0
224 real,
dimension(ngradvar):: max_grad=0.0
225 real,
dimension(ngradvar):: min_grad=0.0
226 character(len=32),
dimension(ngradvar):: grad_var=
"" 228 logical:: force_cubical=.false.
229 logical:: detailed_timer=.false.
232 type(refine_t),
public:: refine
238 SUBROUTINE init (self)
239 class(refine_t):: self
241 integer,
save:: levelmin=0, levelmax=0
242 integer,
save:: ratio=2
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' 252 call trace%begin (
'refine_t%init')
253 call trace%print_id (id)
257 levelmin = shared%levelmin
258 levelmax = shared%levelmax
260 read (io%input, refine_params, iostat=iostat)
261 write (io%output, refine_params)
263 omp_lock%links = .true.
265 omp_lock%links = .false.
273 self%levelmin = levelmin
274 self%levelmax = levelmax
275 shared%levelmin = levelmin
276 shared%levelmax = levelmax
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
295 class(solver_t),
pointer:: patch
296 class(extras_t),
pointer:: extras
298 integer,
save:: itimer=0
302 call trace%begin(
'refine_t%needed', itimer=itimer, detailed_timer=detailed_timer)
306 patch => cast2solver(link%task)
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' 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' 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' 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' 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' 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' 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' 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' 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' 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' 364 call trace%end (itimer, detailed_timer=detailed_timer)
365 END FUNCTION refine_needed
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
380 integer:: refine, id, n_added, l(3), u(3)
381 integer,
save:: itimer=0
383 was_refined = .false.
384 was_derefined = .false.
389 link%task%istep < 3 .or. &
390 link%task%is_set (bits%init_nbors)) &
393 call trace%begin(
'refine_t%check_current', itimer=itimer)
394 patch => cast2solver(link%task)
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
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
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 419 tasklist%verbose = verbose
423 if (.not.
allocated(patch%irefine)) &
424 allocate (patch%irefine(patch%gn(1),patch%gn(2),patch%gn(3)))
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' 439 refine = max(refine,self%refine_needed (tasklist, link))
441 if (verbose > 3)
write (io_unit%output,*) &
442 patch%id, refine, minval(patch%irefine), maxval(patch%irefine),
'mask' 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)
457 patch%check_refine_next = patch%istep + check_interval
458 deallocate (patch%irefine)
460 call trace%end (itimer)
461 END SUBROUTINE check_current
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
472 logical:: was_refined, was_derefined
474 class(solver_t),
pointer:: child
475 class(link_t),
pointer:: nbor, link, nbors, tail
476 class(patch_t),
pointer:: nbpatch
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(:,:,:)
485 call trace%begin (
'refine_t%selective_refine', itimer=itimer, detailed_timer=detailed_timer)
486 was_refined = .false.
487 was_derefined = .false.
489 call link%lock%set (
'selective_refine')
491 tail => tasklist%tail
492 nbor => patch%link%nbor
496 allocate (refined(patch%gn(1),patch%gn(2),patch%gn(3)))
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)
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)))
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
521 call link%lock%unset (
'selective_refine')
525 if (patch%level < self%levelmax .and. refine > 0)
then 526 if (any(patch%irefine > 0))
then 533 size = patch%size/self%ratio
534 pos = patch%position + ([i,j,k]-(self%ratio+1)*0.5)*
size 538 n = patch%mesh%n/self%ratio
539 l = patch%mesh%li + ([i,j,k]-1)*n
544 refine_it = maxval(patch%irefine(l(1):u(1),l(2):u(2),l(3):u(3)))
549 not_already_refined = refined(ii(1),ii(2),ii(3)) == 0
550 if (not_already_refined .and. refine_it > 0)
then 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)
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)
584 if (refine == -1 .and. &
585 patch%level > self%levelmin .and. &
586 patch%check_refine_next > 0 .and. &
587 patch%istep >= 4 .and. &
588 .not.was_refined)
then 591 if (all(patch%irefine(l(1):u(1),l(2):u(2),l(3):u(3)) == -1))
then 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
606 call data_io%update_counters (patch, -1)
607 call patch%log (
'remove')
608 call tasklist%remove_and_reset (link)
609 was_derefined = .true.
613 call trace%end (itimer, detailed_timer=detailed_timer)
614 END SUBROUTINE selective_refine
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
637 class(link_t),
pointer:: old, nbor
638 class(patch_t),
pointer:: old_patch, new_patch, nbor_patch
639 type(patch_t):: patch
641 logical:: overlap, guards
642 integer,
allocatable,
dimension(:,:,:):: trefine
648 if (new%task%level <= self%levelmin+1)
return 650 call trace%begin (
'refine_t%check_support', itimer=itimer, detailed_timer=detailed_timer)
656 new_patch => patch%cast2patch (new%task)
657 if (
associated(new_patch))
then 659 write (stdout,*)
'check_support: id =', &
660 new_patch%id, new_patch%is_set (bits%support)
664 call new_patch%clear (bits%support)
665 call tlist%lock%set (
'check_support')
667 do while (
associated(old))
668 if (old%task%level <= new%task%level-2)
then 669 old_patch => patch%cast2patch (old%task)
671 if (old_patch%get_overlap (new_patch, guards, l, u))
then 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)
679 allocate (trefine(old_patch%gn(1),old_patch%gn(2),old_patch%gn(3)))
681 trefine(l(1):u(1), l(2):u(2), l(3):u(3)) = 1
686 do while (
associated(nbor))
687 if (nbor%task%level == new%task%level-1)
then 688 nbor_patch => patch%cast2patch (nbor%task)
690 if (old_patch%get_overlap (nbor_patch, guards, l, u))
then 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
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
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
725 call tlist%lock%unset (
'check_support')
728 call trace%end (itimer, detailed_timer=detailed_timer)
729 END SUBROUTINE check_support
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
740 class(link_t),
pointer:: nbor1, nbor2, nbor3
741 class(patch_t),
pointer:: task, nbtask
742 type(patch_t):: patch
746 call trace%begin (
'refine_t%need_to_support')
747 task => patch%cast2patch(link%task)
751 call tlist%lock%set (
'need_to_support')
753 do while (
associated(nbor1))
754 if (nbor1%task%level >= link%task%level+2)
then 755 nbtask => patch%cast2patch (nbor1%task)
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
766 call link%lock%set (
'need_to_support')
768 do while (
associated(nbor1))
769 if (nbor1%task%level == link%task%level+1)
then 770 nbtask => patch%cast2patch (nbor1%task)
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
778 call link%lock%unset (
'need_to_support')
779 call tlist%lock%unset (
'need_to_support')
782 refine = maxval(task%irefine(l(1):u(1),l(2):u(2),l(3):u(3)))
784 END FUNCTION need_to_support
792 FUNCTION check_derefine_support (self, tlist, link, check_only)
RESULT (refine)
794 class(refine_t):: self
795 class(list_t):: tlist
796 class(link_t),
pointer:: link
797 logical,
optional:: check_only
800 class(link_t),
pointer:: nbor, nbors
801 integer:: n_nbor1, n_nbor2, n_added
802 integer,
save:: itimer=0
805 if (link%task%level == self%levelmin) &
807 call trace%begin(
'refine_t%check_derefine_support', itimer=itimer, detailed_timer=detailed_timer)
809 call tlist%lock%set (
'check_derefine_support')
811 do while (
associated(nbor))
812 if (nbor%task%level == link%task%level+1)
then 819 write (stdout,*) link%task%id,
'check_derefine_support =', refine
823 call tlist%lock%unset (
'check_derefine_support')
825 call trace%end (itimer, detailed_timer=detailed_timer)
826 END FUNCTION check_derefine_support
832 SUBROUTINE make_child_patch (self, tlist, parent, pos, child_size, &
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
839 class(link_t),
pointer :: child_link, child_nbor
840 integer:: j, count, l(3), u(3), n_added
841 integer,
save:: itimer=0
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')
856 call parent%lock%set (
'parent')
857 allocate (child, source=parent)
858 call parent%lock%unset (
'parent')
859 call child%clone_mem_accounting
861 child%mem_allocated = .false.
865 nullify (child%mem_lock)
867 child%size = child_size
869 child%llc_cart = child%position - &
870 merge(0.5,0.0,child%n>1) * child%size
871 child%llc_nat = child%llc_cart
878 child%n = shared%patch_n
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
892 allocate (child_link)
894 child_link%task => child
895 child_link%parent => parent%link
896 child%link => child_link
897 child%parent => parent
901 child%level = parent%level + 1
902 call child%set (bits%frozen)
904 call child%clear (bits%frozen)
905 child%level = parent%level + 1
906 child%time = parent%time
910 child%iit(child%nt) = child%new
911 child%dt(:) = parent%dt(parent%it) / self%ratio
912 child%t(:) = child%time
913 child%parentid = parent%id
914 child%check_refine_next = 0
915 child%iout = child%time/io%out_time+1
916 child%out_next = child%iout*io%out_time
917 child%dnload_time = -1.0
921 call data_io%update_counters (child, +1)
922 call child%log (
'created')
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
934 if (child%time > io%end_time)
then 935 call child%set (bits%frozen)
945 call child%set (bits%support)
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
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.)
970 call tlist%add_new_link (child_link)
971 call tlist%queue_by_time (child_link)
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))
991 call trace%end (itimer, detailed_timer=detailed_timer)
992 END SUBROUTINE make_child_patch
997 FUNCTION check_cubical (self, patch)
RESULT (refine_it)
998 class(refine_t):: self
999 class(solver_t),
pointer:: patch
1002 if (all(patch%n==shared%patch_n))
then 1008 END FUNCTION check_cubical
1013 FUNCTION check_jeans (self, patch)
RESULT(refine_it)
1016 class(refine_t):: self
1017 class(solver_t),
pointer:: patch
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(:,:,:)
1027 call trace%begin (
'refine_t%check_jeans')
1028 call self%sanity_check (min_jeans, max_jeans,
'check_Jeans')
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)))
1035 pg = patch%gas_pressure() + tiny(1.0)
1037 call mpi%abort(
"Unrecognised patch type. Abort!")
1039 dmin = patch%fminval(d)
1040 if (dmin <= 0.0)
then 1042 write (io%output,*) &
1043 'check_Jeans ERROR: id, dmin, imin =', patch%id, dmin, imin
1044 jeans_num = min_jeans
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))
1049 test = sqrt(pg)/(c*d)
1055 jeans_num = minval(test, mask=patch%irefine <= 0)
1062 where (test < min_jeans)
1067 else where (test < max_jeans)
1068 patch%irefine = max(patch%irefine, 0)
1072 if (jeans_num > max_jeans) refine_it = -1
1073 if (jeans_num < min_jeans) refine_it = +1
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), &
1088 END FUNCTION check_jeans
1093 FUNCTION check_compress (self, patch)
RESULT(refine_it)
1096 class(refine_t):: self
1097 class(solver_t),
pointer:: patch
1100 integer:: imin(3), l(3), u(3)
1102 real(kind=KindScalarVar),
dimension(:,:,:),
allocatable:: test
1104 integer,
save:: itimer=0
1106 call trace%begin (
'refine_t%check_compress', itimer=itimer, detailed_timer=detailed_timer)
1109 allocate (test(patch%gn(1),patch%gn(2),patch%gn(3)))
1112 call patch%compression_magnitude (test)
1114 call mpi%abort(
"Unrecognised patch type. Abort!")
1116 test = minval(patch%mesh%d)*test
1117 w_num = maxval(test, mask=patch%irefine <= 0)
1124 where (test > max_compress)
1129 else where (test > min_compress)
1130 patch%irefine = max(patch%irefine, 0)
1133 if (w_num > max_compress) refine_it = +1
1134 if (w_num < min_compress) refine_it = -1
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), &
1146 call trace%end (itimer, detailed_timer=detailed_timer)
1147 END FUNCTION check_compress
1152 FUNCTION check_vorticity (self, patch)
RESULT(refine_it)
1155 class(refine_t):: self
1156 class(solver_t),
pointer:: patch
1159 integer:: imin(3), l(3), u(3)
1161 real(kind=KindScalarVar),
dimension(:,:,:),
allocatable:: test
1163 integer,
save:: itimer=0
1165 call trace%begin (
'refine_t%check_vorticity', itimer=itimer, detailed_timer=detailed_timer)
1168 allocate (test(patch%gn(1),patch%gn(2),patch%gn(3)))
1171 call patch%vorticity_magnitude (test)
1173 call mpi%abort(
"Unrecognised patch type. Abort!")
1175 test = minval(patch%mesh%d)*test
1176 w_num = maxval(test, mask=patch%irefine <= 0)
1183 where (test > max_vorticity)
1188 else where (test > min_vorticity)
1189 patch%irefine = max(patch%irefine, 0)
1192 if (w_num > max_vorticity) refine_it = +1
1193 if (w_num < min_vorticity) refine_it = -1
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), &
1205 call trace%end (itimer, detailed_timer=detailed_timer)
1206 END FUNCTION check_vorticity
1211 FUNCTION shock_and_cd_detector (self, patch)
RESULT (refine)
1213 class(refine_t):: self
1214 class(solver_t),
pointer:: patch
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)
1226 call trace%begin(
'refine_t%shocks_and_cd_detector')
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))
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()
1243 call mpi%abort(
"Unrecognised patch type. Abort!")
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 )
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)
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
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
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 )
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
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
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
1305 nullify (mx, my, mz, d)
1309 END FUNCTION shock_and_cd_detector
1314 FUNCTION gradient_detector (self, patch)
RESULT (refine)
1317 class(refine_t):: self
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']
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)))
1338 if (mx%n > 1) ione = 1
1339 if (my%n > 1) jone = 1
1340 if (mz%n > 1) kone = 1
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)
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)
1359 call mpi%abort(
"Unrecognised patch type. Abort!")
1361 is_unsigned = .false.
1362 do i=1,
size(unsigned)
1363 if (grad_var(n) == unsigned(i))
then 1364 is_unsigned = .true.
1375 if (is_unsigned)
then 1379 dq1(i,j,k) = abs( q(ip1,j,k) - q(im1,j,k) ) &
1380 / ( max( q(ip1,j,k), q(im1,j,k) ) &
1382 dq2(i,j,k) = abs( q(i,jp1,k) - q(i,jm1,k) ) &
1383 / ( max( q(i,jp1,k), q(i,jm1,k) ) &
1385 dq3(i,j,k) = abs( q(i,j,kp1) - q(i,j,km1) ) &
1386 / ( max( q(i,j,kp1), q(i,j,km1) ) &
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) )
1402 do k=mz%li+kone,mz%ui-kone
1405 do j=my%li+jone,my%ui-jone
1408 do i=mx%li+ione,mx%ui-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
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
1427 deallocate (dq1, dq2, dq3, q)
1429 END FUNCTION gradient_detector
1436 FUNCTION refine_region (self, patch)
RESULT (refine)
1438 class(refine_t) :: self
1439 class(solver_t),
pointer :: patch
1441 class(
mesh_t),
pointer:: m1, m2, m3, m(:)
1442 logical :: intersect
1443 logical :: refine_it
1444 integer :: ix, iy, iz
1446 real,
dimension(3) :: p_llc, p_urc, b_llc, b_urc
1448 call trace%begin(
'refine_t%refine_on_box')
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
1463 intersect = all((p_urc - b_llc) > 0) .and. all((b_urc - p_llc) > 0)
1464 if (.not. intersect)
then 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) )
1480 refine = merge(1,0,refine_it)
1482 END FUNCTION refine_region
1487 SUBROUTINE sanity_check (self, min, max, label, reverse)
1488 class(refine_t):: self
1490 character(len=*):: label
1491 logical,
optional:: reverse
1493 if (max == 0.0)
then 1494 if (
present(reverse))
then 1495 min = max/self%safe_ratio
1497 max = min*self%safe_ratio
1499 else if (max < min*2.0)
then 1500 if (
present(reverse))
then 1506 write (stdout,*)
'WARNING: '//trim(label)//
' max value set to', max
1509 END SUBROUTINE sanity_check
1514 FUNCTION cast2solver (task)
RESULT (solver)
1515 class(task_t),
pointer:: task
1516 class(solver_t),
pointer:: solver
1523 END FUNCTION cast2solver
Each thread uses a private timer data type, with arrays for start time and total time for each regist...
download_link: takes care of downloads to linktask same: called for patches on the same level differ:...
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.
Define code units, in terms of (here) CGS units.
Interface from gpatch_mod to a choice of binary data I/O methods, controlled by the iomethod text str...
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
This module handles checking max change between neighboring points. Each instance of it needs an inde...
This module is a placeholder for shared information, solving the problem of making a piece of informa...
Fundamental constants in CGS and SI units.
Template module for mesh.
Module with list handling for generic class task_t objects.
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...
Template module for tasks.