7 USE iso_fortran_env
, only: int8
28 integer,
parameter:: ndim = 3
30 real(kind=KindScalarVar),
dimension(:,:,:,:,:,:),
pointer:: mem => null()
31 real(kind=KindScalarVar),
dimension(:,:,:,:,:),
pointer:: vgas => null()
32 real(kind=KindScalarVar),
dimension(:,:,:,:),
pointer:: density => null()
33 real(kind=KindScalarVar),
dimension(:,:,:,:),
pointer:: pressure => null()
39 real(kind=KindScalarVar),
dimension(:,:,:,:),
allocatable:: force_per_unit_mass
40 real(kind=KindScalarVar),
dimension(:,:,:,:),
allocatable:: force_per_unit_volume
41 real(kind=KindScalarVar),
dimension(:,:,:),
allocatable:: heating_per_unit_mass
42 real(kind=KindScalarVar),
dimension(:,:,:),
allocatable:: heating_per_unit_volume
43 integer(kind=int8) ,
dimension(:,:,:),
allocatable:: filled
49 class(
mesh_t),
dimension(:),
pointer:: mesh => null()
51 real(8):: pdf_next=0.0
60 integer:: check_refine_next=0
62 integer:: not_filled=0
63 integer:: refine_ratio=2
65 logical:: mem_allocated=.false.
66 logical:: do_pebbles=.false.
67 logical:: staggered=.true.
68 logical:: no_mans_land=.false.
69 logical:: use_data_hub=.false.
70 logical:: manual_refine=.false.
71 logical:: get_umax_location=.false.
72 logical,
dimension(:),
pointer:: unsigned => null()
73 logical,
dimension(:),
pointer:: pervolume => null()
74 class(
link_t),
pointer:: link
76 character(len=32):: eos=
'ideal' 77 character(len=32):: opacity=
'none' 78 integer:: file_offset(4)
80 real:: safety_factor=1.0
82 logical,
allocatable:: mask(:,:,:)
84 integer(int8),
allocatable:: irefine(:,:,:)
85 type(
lock_t),
pointer:: mem_lock(:) => null()
93 integer:: n(3)=0, ng(3)=0
94 integer:: li(3), ui(3), gn(3)
98 real(8):: llc_cart(3)=0d0
99 real(8):: llc_nat(3)=0d0
100 real(8):: centre_cart(3)=0d0
101 real(8):: centre_nat(3)=0d0
102 integer :: mesh_type = 1
103 real(8):: grid_vel(3) = 0d0
107 procedure,
nopass:: cast2patch
109 procedure:: init_default
110 procedure:: init_bdries
113 procedure:: clone_mem_accounting
115 procedure:: allocate_mesg
117 procedure,
nopass:: unpack_id
120 procedure,
private:: contains1
121 procedure,
private:: contains2
122 generic::
contains => contains1, contains2
128 procedure:: myposition
129 procedure:: myposition_staggered
130 procedure:: periodic_grid
131 procedure:: extrapolate_guards
132 procedure:: is_periodic
133 procedure,
private :: make_periodic4
134 procedure,
private :: make_periodic8
135 generic,
public :: make_periodic => make_periodic4 , make_periodic8
136 procedure,
private :: fsum4
137 procedure,
private :: fsum8
138 generic,
public :: fsum => fsum4, fsum8
139 procedure,
private :: faver4
140 procedure,
private :: faver8
141 generic,
public :: faver => faver4, faver8
142 procedure,
private :: fminval4
143 procedure,
private :: fminval8
144 procedure,
private :: fminvali
145 generic,
public :: fminval => fminval4, fminval8, fminvali
146 procedure,
private :: fmaxval4
147 procedure,
private :: fmaxval8
148 procedure,
private :: fmaxvali
149 generic,
public :: fmaxval => fmaxval4, fmaxval8, fmaxvali
150 generic,
public :: fminvalvar => fminval4, fminval8
151 generic,
public :: fmaxvalvar => fmaxval4, fmaxval8
152 procedure:: courant_condition
153 procedure:: check_density
154 procedure:: check_nan
155 procedure:: custom_refine
156 procedure:: patch_to_header
157 procedure:: header_to_patch
158 procedure,
private:: stats_patch
159 procedure,
private:: stats_scalar
160 procedure,
private:: stats_vector
161 generic,
public:: stats => stats_patch, stats_scalar, stats_vector
162 procedure:: mapvectors
163 procedure:: beyond_patch_edge
164 procedure:: myfloatindex_scalar
165 procedure:: log_interpolate
166 procedure:: time_interval
167 procedure:: distance => distance_curvilinear
168 procedure:: index_stagger
169 procedure:: interpolate
170 procedure:: interpolate_specific
171 procedure:: index_space
172 procedure:: index_space_of
173 procedure:: index_space2
174 procedure:: index_only
175 procedure:: index_only2
176 procedure:: count_edges
177 procedure:: update_position
178 procedure:: get_overlap
179 procedure:: init_level
184 integer,
public,
save:: n_header = (14*4 + 4*8 + 5*8*3 + 1*8 + 2*8*8 + 1*4*8 + 2*4*3)/4
205 real(8):: position(3)
206 real(8):: velocity(3)
211 real(8),
dimension(8):: t
212 real(8),
dimension(8):: dt
213 integer,
dimension(8):: iit
225 integer,
save:: print_every=0, iprint=0
226 logical,
save:: extrapolate_gz=.false.
227 logical,
save:: do_check_nan=.false.
228 logical,
save:: zero_order_extrap=.true.
235 SUBROUTINE patch_to_header (self, head)
240 head%n_header = n_header
241 head%n_data = product(self%gn)*self%nv
243 head%status = self%status
244 head%rank = self%rank
245 head%istep = self%istep
246 head%iout = self%iout
247 head%level = self%level
248 head%time = self%time
249 head%out_next = self%out_next
250 head%print_next = self%print_next
251 head%dtime = self%dtime
252 head%position = self%position
253 head%velocity = self%velocity
254 head%size = self%size
258 head%t(1:nt) = self%t
259 head%dt(1:nt) = self%dt
260 head%iit(1:nt) = self%iit
262 head%send_time = wallclock()
263 if (io%verbose > 1) &
264 write (io_unit%log,
'(f12.6,3x,a,i7,i4,l3,2i4,f8.4)') &
265 wallclock(),
'patch_to_header: id, nq =', self%id, &
266 self%nq, self%is_set(bits%swap_request), &
267 mpi_mesg%sent_list%n, mpi_mesg%recv_list%n
268 END SUBROUTINE patch_to_header
273 SUBROUTINE header_to_patch (self, head)
283 if (iand(head%status, bits%boundary) /= 0)
then 284 head%status = ior(head%status, bits%virtual)
285 head%status = iand(head%status, not(bits%boundary))
287 write (io_unit%mpi, *) head%id,
'swapping to virtual' 288 else if (iand(head%status, bits%virtual) /= 0)
then 289 head%status = ior(head%status, bits%boundary)
290 head%status = iand(head%status, not(bits%virtual))
292 write (io_unit%mpi, *) head%id,
'swapping to boundary' 294 head%status = iand(head%status, not(bits%ready+bits%busy))
295 if (head%id /= self%id)
then 296 call io%abort(
'wrong message ID unpacked')
301 if ((self%istep > 0) .and. (iand(head%status, bits%remove) == 0))
then 302 if (head%istep > self%istep+1)
then 303 write(io_unit%log,1) &
304 'WARNING: early arrival in time reversal id, rank =', &
305 self%id, mpi%rank, head%istep, self%istep, head%time, self%time
306 1
format(a,4i6,2f12.6)
308 else if (head%istep < self%istep+1)
then 309 write(io_unit%log,1) &
310 'WARNING: late arrival in time reversal id, rank =', &
311 self%id, mpi%rank, head%istep, self%istep, head%time, self%time
316 self%status = head%status
317 self%rank = head%rank
318 self%istep = head%istep
319 self%iout = head%iout
320 self%level = head%level
321 self%time = head%time
322 self%out_next = head%out_next
323 self%print_next = head%print_next
324 self%dtime = head%dtime
325 self%position = head%position
326 self%velocity = head%velocity
327 self%size = head%size
330 new = mod(head%it,head%nt) + 1
336 self%t = head%t(1:nt)
337 self%dt = head%dt(1:nt)
338 self%iit = head%iit(1:nt)
341 write (io_unit%log,*) omp%thread,
'header_to_patch WARNING: nt non-positive', nt
345 self%latency = wallclock() - head%send_time
347 timer%latency%max = max(timer%latency%max,self%latency)
349 timer%latency%aver = timer%latency%aver+self%latency
351 timer%latency%n = timer%latency%n+1
353 END SUBROUTINE header_to_patch
358 FUNCTION cast2patch (task)
RESULT(patch)
359 class(
task_t),
pointer:: task
360 class(
patch_t),
pointer:: patch
367 call io%abort (
'patch_t%cast: failed to cast a task to patch_t')
369 END FUNCTION cast2patch
374 FUNCTION task2patch (task)
RESULT (patch)
376 class(
patch_t),
pointer:: patch
378 patch => cast2patch(task)
379 END FUNCTION task2patch
386 SUBROUTINE init (self)
388 character(len=120):: id = &
389 '$Id: f01f80ce28846df45c29d0d9c93488206e197f17 $ tasks/patch_mod.f90' 391 call trace%begin (
'patch_t%init')
392 call trace%print_id (id)
421 SUBROUTINE init_default (self)
424 integer,
save:: n(3)=16, ng(3)=2, nv=5, nt=5, nw, max_files=1000
425 logical,
save:: no_mans_land=.false., use_data_hub=.false., limit_dtime=.false.
426 real,
save:: u_fixed=-1.0, grace=0.0
427 real(8),
save:: end_time=1.0, out_time=0.1, print_time=0.1, dt_fixed=-1.0
428 real(8):: position(3), size(3)
429 integer:: iostat, i, n_alloc=0, n_ext=0
430 namelist /patch_params/ n, nv, ng, nt, u_fixed, dt_fixed, &
431 grace, extrapolate_gz, no_mans_land, use_data_hub, do_check_nan, &
432 zero_order_extrap, limit_dtime
433 namelist /out_params/ out_time, end_time, print_time, print_every, &
435 logical,
save:: first_time = .true.
437 call trace%begin (
'patch_t%init_default')
438 call self%task_t%init
446 if (any(self%n /= 0)) n = self%n
447 if (any(self%ng /= 0)) ng = self%ng
448 if ( self%nv /= 0 ) nv = self%nv
449 if ( self%nt /= 0 ) nt = self%nt
452 read (io%input, patch_params)
453 if (dt_fixed > 0.0) grace=max(grace,1e-6)
454 if (io%master)
write (*, patch_params)
456 read (io%input, out_params)
457 if (io%master)
write (*, out_params)
466 if (any(self%n == 0)) self%n = n
467 if (any(self%ng == 0)) self%ng = ng
468 if ( self%nv == 0 ) self%nv = nv
469 if ( self%nt == 0 ) self%nt = nt
473 io%end_time = end_time
474 io%out_time = out_time
475 io%print_time = print_time
479 self%restart = io%restart
481 self%dt_fixed = dt_fixed
482 self%u_fixed = u_fixed
483 self%max_files = max_files
484 self%use_data_hub = use_data_hub
485 self%no_mans_land = no_mans_land
486 self%print_next = print_time
487 self%limit_dtime = limit_dtime
489 END SUBROUTINE init_default
491 SUBROUTINE init_bdries (self)
493 real(8):: position(3), limit(3)
497 call trace%begin (
'patch_t%init_bdries')
498 position = self%position - self%origin - 0.5d0 * self%box
499 limit = 0.5d0 * (self%box - self%size - self%ds)
500 if (.not.self%periodic(1))
then 501 if (position(1) <= -limit(1))
then 502 call self%boundaries%set (bits%xl)
503 self%mesh(1)%lower_boundary = .true.
505 if (position(1) >= +limit(1))
then 506 call self%boundaries%set (bits%xu)
507 self%mesh(1)%upper_boundary = .true.
510 if (.not.self%periodic(2))
then 511 if (position(2) <= -limit(2))
then 512 call self%boundaries%set (bits%yl)
513 self%mesh(2)%lower_boundary = .true.
515 if (position(2) >= +limit(2))
then 516 call self%boundaries%set (bits%yu)
517 self%mesh(2)%upper_boundary = .true.
520 if (.not.self%periodic(3))
then 521 if (position(3) <= -limit(3))
then 522 call self%boundaries%set (bits%zl)
523 self%mesh(3)%lower_boundary = .true.
525 if (position(3) >= +limit(3))
then 526 call self%boundaries%set (bits%zu)
527 self%mesh(3)%upper_boundary = .true.
531 END SUBROUTINE init_bdries
537 SUBROUTINE setup (self, position, size, n, ng, nt, nv, nw)
539 real(8),
optional :: position(3), size(3)
540 integer,
optional :: n(3), ng(3)
541 integer,
optional :: nt, nv, nw
544 character(len=4) kind
547 call trace%begin (
'patch_t%setup')
548 call self%init_default
549 if (
present(position)) self%position = position
550 if (
present(
size )) self%size =
size 551 if (
present(n )) self%n = n
552 if (
present(ng )) self%ng = ng
553 if (
present(nt )) self%nt = nt
554 if (
present(nv )) self%nv = nv
555 if (
present(nw )) self%nw = nw
562 shared%patch_n = self%n
564 self%llc_cart = self%position - 0.5 * self%size
565 self%llc_nat = self%llc_cart
569 if (all(self%size==self%box))
then 571 write (io_unit%log,*) mpi%rank,
'task', self%id,
' is root grid' 572 call self%set(bits%root_grid)
582 self%llc_cart = self%position - 0.5_8*self%size
583 if (self%mesh_type == mesh_types%Cartesian)
then 584 self%llc_nat = self%llc_cart
585 self%centre_nat = self%position
591 call omp%set_stacksize (self%nv)
592 call meshfactory(self%mesh_type, self%mesh, n=self%n, ng=self%ng, s=self%size, &
593 p=self%position, llc_native=self%llc_nat, nml=self%no_mans_land)
597 allocate(self%mesh(i)%h(self%nv))
599 if (self%kind ==
'zeus_mhd_patch' .and. &
600 self%staggered .and. self%n(i) /= 1) self%mesh(i)%gn = self%mesh(i)%gn + 1
601 if (self%kind ==
'zeus_tw_mhd_patch' .and. &
602 self%staggered .and. self%n(i) /= 1) self%mesh(i)%gn = self%mesh(i)%gn + 1
604 self%ds = self%mesh%d
605 self%ncell = self%mesh%nc
606 self%mesh%id = self%id
607 self%mesh%b = self%box
608 self%mesh%llc_cart = self%llc_cart
609 self%centre_nat = [self%mesh(1)%centre_nat,self%mesh(2)%centre_nat,self%mesh(3)%centre_nat]
610 self%centre_cart = currenttocartesian(self%mesh_type, &
611 self%centre_nat(1), self%centre_nat(2), self%centre_nat(3))
615 self%gn = self%mesh%gn
616 self%li = self%mesh%li
617 self%ui = self%mesh%ui
618 self%gsize = (self%size*self%gn/self%n)
619 self%offset = (self%mesh%lb+self%mesh%ub)/2.0
620 if (self%id == io%id_debug)
then 621 print 1,
'patch_t%init: position =', self%position
622 print 1,
'patch_t%init: m%p =', self%mesh%p
623 print 1,
'patch_t%init: offset =', self%offset
624 print 1,
'patch_t%init: m%o =', self%mesh%o
625 print 1,
'patch_t%init: centre_nat =', self%centre_nat
626 print 1,
'patch_t%init: centre_cart =', self%centre_cart
627 print 1,
'patch_t%init: llc_nat =', self%llc_nat
628 print 1,
'patch_t%init: m%llc_nat =', self%mesh%llc_nat
629 print 1,
'patch_t%init: llc_cart =', self%llc_cart
630 print 1,
'patch_t%init: m%llc_cart =', self%mesh%llc_cart
631 1
format(a,1p,3e15.6)
638 if (io%verbose>2)
then 640 write (io_unit%log,
'(1x,a,i7,2(1p,3e12.4,2x))')
'patch_mod::init: id, size, pos =', &
641 self%id, self%size, self%position
643 write (io_unit%log,*)
'id, LEVEL: ', self%id, self%level, maxval(self%box), minval(self%ds)
649 if (.not.self%mem_allocated)
then 650 allocate (self%t(self%nt), self%dt(self%nt), self%iit(self%nt))
653 allocate (self%mem(self%mesh(1)%gn, &
656 self%nv,self%nt,self%nw))
657 allocate (self%mem_lock(self%nt))
659 write (kind,
'("mem",i1)') i
660 call self%mem_lock(i)%init (kind)
662 self%mem_allocated = .true.
663 call io%bits_mem (storage_size(self%mem), product(shape(self%mem)),
'mem')
664 allocate (self%unsigned (self%nv))
665 allocate (self%pervolume(self%nv))
666 self%unsigned (:) = .false.
667 self%pervolume(:) = .false.
674 self%new = min(2,self%nt)
675 self%iit(self%nt) = self%new
676 self%new = self%iit(self%nt)
677 self%t(self%it) = 0.0
678 self%wallclock = wallclock()
679 n_header = storage_size(header)/32
680 if (n_header*4 /= storage_size(header)/8)
then 681 write(io%output,*) n_header*4, storage_size(header)/8
682 call mpi%abort (
'n_header is incorrect')
684 if (all(self%size == self%box))
then 685 call self%set (bits%root_grid)
686 if (io%verbose > 1)
write(io%output,*)
'set bits%root_grid id =', self%id
694 SUBROUTINE dealloc (self)
698 call trace%begin (
'patch_t%dealloc')
699 if (
associated(self%mesh)) &
700 call meshrecycler (self%mesh)
701 call io%bits_mem (-storage_size(self%mem), product(shape(self%mem)),
'-mem')
702 if (
associated(self%mem))
deallocate (self%mem)
704 self%mem_allocated = .false.
705 if (
allocated(self%force_per_unit_mass))
then 706 call io%bits_mem (-storage_size(self%force_per_unit_mass), &
707 product(shape(self%force_per_unit_mass)),
'-fpm')
708 deallocate (self%force_per_unit_mass)
710 if (
allocated(self%force_per_unit_volume))
then 711 call io%bits_mem (-storage_size(self%force_per_unit_volume), &
712 product(shape(self%force_per_unit_volume)),
'-fpv')
713 deallocate (self%force_per_unit_volume)
715 if (
allocated(self%heating_per_unit_mass))
then 716 call io%bits_mem (-storage_size(self%heating_per_unit_mass), &
717 product(shape(self%heating_per_unit_mass)),
'-qpm')
718 deallocate (self%heating_per_unit_mass)
720 if (
allocated(self%heating_per_unit_volume))
then 721 call io%bits_mem (-storage_size(self%heating_per_unit_volume), &
722 product(shape(self%heating_per_unit_volume)),
'-qpv')
723 deallocate (self%heating_per_unit_volume)
725 if (
associated(self%vgas))
then 726 call io%bits_mem (-storage_size(self%vgas), &
727 product(shape(self%vgas)),
'-vgas')
728 deallocate (self%vgas)
730 if (
associated(self%density))
then 731 call io%bits_mem (-storage_size(self%density), &
732 product(shape(self%density)),
'-density')
733 deallocate (self%density)
735 if (
associated(self%pressure))
then 736 call io%bits_mem (-storage_size(self%pressure), &
737 product(shape(self%pressure)),
'-pressure')
738 deallocate (self%pressure)
740 if (
allocated(self%mask))
then 741 call io%bits_mem (-storage_size(self%mask), &
742 product(shape(self%mask)),
'-mask')
743 deallocate (self%mask)
745 if (
allocated(self%irefine))
then 746 call io%bits_mem (-storage_size(self%irefine), &
747 product(shape(self%irefine)),
'-irefine')
748 deallocate (self%irefine)
750 call self%task_t%dealloc
752 END SUBROUTINE dealloc
754 SUBROUTINE clone_mem_accounting (self)
756 if (
allocated(self%mask)) &
757 call io%bits_mem (storage_size(self%mask), &
758 product(shape(self%mask)),
'mask')
759 if (
allocated(self%irefine)) &
760 call io%bits_mem (storage_size(self%irefine), &
761 product(shape(self%irefine)),
'irefine')
762 if (
allocated(self%force_per_unit_mass)) &
763 call io%bits_mem (storage_size(self%force_per_unit_mass), &
764 product(shape(self%force_per_unit_mass)),
'fpm')
765 if (
allocated(self%force_per_unit_volume)) &
766 call io%bits_mem (storage_size(self%force_per_unit_volume), &
767 product(shape(self%force_per_unit_volume)),
'fpv')
768 if (
allocated(self%heating_per_unit_mass)) &
769 call io%bits_mem (storage_size(self%heating_per_unit_mass), &
770 product(shape(self%heating_per_unit_mass)),
'qpm')
771 if (
allocated(self%heating_per_unit_volume)) &
772 call io%bits_mem (storage_size(self%heating_per_unit_volume), &
773 product(shape(self%heating_per_unit_volume)),
'qpv')
774 END SUBROUTINE clone_mem_accounting
782 SUBROUTINE rotate (self)
784 integer:: iv, nv, new
786 integer,
save:: itimer=0
790 if (self%rotated)
return 791 call trace%begin (
'patch_t%rotate', itimer=itimer)
792 if (self%is_periodic())
then 793 call self%periodic_grid
794 else if (extrapolate_gz)
then 795 call self%extrapolate_guards
800 if (io%verbose>1)
then 801 write (io_unit%log,
'(a,i6,2i3,2e13.6,2x,a,1p,6e10.2)')
'update: id,j,k,l,n,dt,t', &
802 self%id, self%it, self%new, self%dtime, self%time,
'minmax:', &
803 minval(self%mem(:,:,:,1,self%it ,self%nw)), &
804 maxval(self%mem(:,:,:,1,self%it ,self%nw)), &
805 minval(self%mem(:,:,:,1,self%new,1)), &
806 maxval(self%mem(:,:,:,1,self%new,1))
812 call self%update_position
816 call self%task_t%rotate
817 call trace%end (itimer)
818 END SUBROUTINE rotate
824 FUNCTION myposition (self, ii, iv, pp)
RESULT(out)
826 integer,
intent(in) :: ii(3)
827 integer,
optional :: iv
828 real,
optional :: pp(3)
830 real(8) :: out(3), cellpos(3), h(3)
832 if (
present(pp))
then 833 cellpos = self%position + (
real(ii)-self%offset + pp)*self%ds
835 cellpos = self%position + (
real(ii)-self%offset)*self%ds
837 if (
present(iv))
then 838 h = self%index_stagger(iv)
839 cellpos = cellpos + h*self%mesh%d
841 if (self%id==io%id_debug .and. io%verbose>1) &
842 write(io%output,
'(a,3i4,4(3x,3f10.2))') &
843 'mypos:',ii,self%position,self%offset,self%ds,cellpos
845 if (self%mesh_type == mesh_types%cylindrical)
then 847 else if (self%mesh_type == mesh_types%spherical)
then 848 cellpos = sphericaltocartesian(cellpos(1), cellpos(2), cellpos(3))
851 END FUNCTION myposition
857 FUNCTION myposition_staggered (self, ii, idir, h)
RESULT(out)
859 integer,
intent(in):: ii(3), idir
860 real(8):: out, cellpos(3)
861 real(8),
optional :: h
870 cellpos = (
real(ii-self%li,kind=8) + merge(0.5 + hlocal,0.0,self%n(idir)>1)) &
871 * self%ds + self%llc_nat
873 if (self%mesh_type == mesh_types%cylindrical)
then 875 else if (self%mesh_type == mesh_types%spherical)
then 876 cellpos = sphericaltocartesian(cellpos(1), cellpos(2), cellpos(3))
880 END FUNCTION myposition_staggered
883 FUNCTION myfloatindex_scalar (self, cpos, idir, h)
RESULT(out)
885 real(8),
intent(in) :: cpos(3)
886 integer,
intent(in) :: idir
887 real(8),
intent(in) :: h
888 real(8) :: pp(3), ppCart(3), out, intermezzo
891 if (self%mesh_type == mesh_types%Cartesian)
then 894 else if (self%mesh_type == mesh_types%spherical)
then 895 ppcart = cpos - self%position
896 pp = cartesiantospherical(ppcart(1), ppcart(2), ppcart(3))
897 else if (self%mesh_type == mesh_types%cylindrical)
then 898 ppcart = cpos - self%position
899 pp = cartesiantocylindrical(ppcart(1), ppcart(2), ppcart(3))
902 intermezzo = (pp(idir) - self%llc_nat(idir)) / self%ds(idir) + self%li(idir) &
903 - merge(0.5d0 + h,0.0d0,self%n(idir)>1)
906 END FUNCTION myfloatindex_scalar
912 SUBROUTINE periodic_grid (self, only)
914 integer,
optional:: only
916 integer:: ix, iy, iz, i(3), j(3), lb(3), ub(3), li(3), ui(3), n(3)
917 logical:: periodic(3)
919 call trace%begin (
'patch_t%periodic_grid')
920 periodic = self%periodic .and. self%size+self%ds > self%box
921 if (.not.any(periodic))
call mpi%abort (
'patch_t%periodic_grid')
926 where (self%mesh%lower_boundary)
929 where (self%mesh%upper_boundary)
937 if (all(i >= li .and. i <= ui)) cycle
939 j = modulo(i-li,n)+li
946 where (.not.periodic)
955 where (self%mesh%lower_boundary .and. i < self%mesh%li+4)
958 where (self%mesh%upper_boundary .and. i > self%mesh%ui-4)
961 if (
present(only))
then 962 self%mem(ix,iy,iz,only,self%it ,1) = self%mem(j(1),j(2),j(3),only,self%it ,1)
963 self%mem(ix,iy,iz,only,self%new,1) = self%mem(j(1),j(2),j(3),only,self%new,1)
965 self%mem(ix,iy,iz,:,self%it ,1) = self%mem(j(1),j(2),j(3),:,self%it ,1)
966 self%mem(ix,iy,iz,:,self%new,1) = self%mem(j(1),j(2),j(3),:,self%new,1)
971 if (io%verbose > 1)
write(io%output,*) &
972 'periodic_grid: id, it, new =', self%id, self%it, self%new
974 END SUBROUTINE periodic_grid
977 FUNCTION is_periodic(self)
RESULT(out)
981 out = any(self%periodic .and. self%size+self%ds > self%box)
982 END FUNCTION is_periodic
987 SUBROUTINE make_periodic4 (patch, f)
989 class(
mesh_t),
pointer:: m(:)
991 integer:: ix, iy, iz, n(3)
992 logical:: periodic(3)
994 periodic = patch%periodic .and. patch%size+patch%ds > patch%box
1000 if (n(1) > 1 .and. periodic(1))
then 1001 do iz=m(3)%li,m(3)%ui
1002 do ix=m(1)%lb,m(1)%lo
1003 do iy=m(2)%li,m(2)%ui
1004 f(ix,iy,iz) = f(ix+n(1),iy,iz)
1007 do ix=m(1)%uo,m(1)%ub
1008 do iy=m(2)%li,m(2)%ui
1009 f(ix,iy,iz) = f(ix-n(1),iy,iz)
1017 if (n(2) > 1 .and. periodic(2))
then 1018 do iz=m(3)%li,m(3)%ui
1019 do iy=m(2)%lb,m(2)%lo
1020 do ix=m(1)%lb,m(1)%ub
1021 f(ix,iy,iz) = f(ix,iy+n(2),iz)
1024 do iy=m(2)%uo,m(2)%ub
1025 do ix=m(1)%lb,m(1)%ub
1026 f(ix,iy,iz) = f(ix,iy-n(2),iz)
1034 if (n(3) > 1 .and. periodic(3))
then 1035 do iy=m(2)%lb,m(2)%ub
1036 do iz=m(3)%lb,m(3)%lo
1037 do ix=m(1)%lb,m(1)%ub
1038 f(ix,iy,iz) = f(ix,iy,iz+n(3))
1041 do iz=m(3)%uo,m(3)%ub
1042 do ix=m(1)%lb,m(1)%ub
1043 f(ix,iy,iz) = f(ix,iy,iz-n(3))
1048 END SUBROUTINE make_periodic4
1053 SUBROUTINE make_periodic8 (patch, f)
1055 class(
mesh_t),
pointer:: m(:)
1057 integer:: ix, iy, iz, n(3)
1058 logical:: periodic(3)
1060 periodic = patch%periodic .and. patch%size+patch%ds > patch%box
1066 if (n(1) > 1 .and. periodic(1))
then 1067 do iz=m(3)%li,m(3)%ui
1068 do ix=m(1)%lb,m(1)%lo
1069 do iy=m(2)%li,m(2)%ui
1070 f(ix,iy,iz) = f(ix+n(1),iy,iz)
1073 do ix=m(1)%uo,m(1)%ub
1074 do iy=m(2)%li,m(2)%ui
1075 f(ix,iy,iz) = f(ix-n(1),iy,iz)
1083 if (n(2) > 1 .and. periodic(2))
then 1084 do iz=m(3)%li,m(3)%ui
1085 do iy=m(2)%lb,m(2)%lo
1086 do ix=m(1)%lb,m(1)%ub
1087 f(ix,iy,iz) = f(ix,iy+n(2),iz)
1090 do iy=m(2)%uo,m(2)%ub
1091 do ix=m(1)%lb,m(1)%ub
1092 f(ix,iy,iz) = f(ix,iy-n(2),iz)
1100 if (n(3) > 1 .and. periodic(3))
then 1101 do iy=m(2)%lb,m(2)%ub
1102 do iz=m(3)%lb,m(3)%lo
1103 do ix=m(1)%lb,m(1)%ub
1104 f(ix,iy,iz) = f(ix,iy,iz+n(3))
1107 do iz=m(3)%uo,m(3)%ub
1108 do ix=m(1)%lb,m(1)%ub
1109 f(ix,iy,iz) = f(ix,iy,iz-n(3))
1114 END SUBROUTINE make_periodic8
1121 SUBROUTINE extrapolate_guards (self)
1124 integer :: it0, it1, it2, nt, iit(self%nt)
1125 integer,
dimension(3) :: lb, lo, li, ub, uo, ui, i
1127 integer,
save :: itimer=0
1128 real(8) :: t(self%nt)
1130 call trace%begin (
'patch_t%extrapolate_guards', 3, itimer=itimer)
1132 call self%timeslots (iit, t)
1145 where (self%mesh%lower_boundary)
1148 where (self%mesh%upper_boundary)
1151 if (io%id_debug>0.and.(self%id==io%id_debug.or.abs(self%id-io%id_debug)==27))
then 1152 write(io%output,
'(i5,2x,a,6i7)') self%id, &
1153 'MK extrapolating guard zone values', it2, it1, it0, self%it, self%new
1156 if (self%unsigned(iv))
then 1160 if (self%gn(3) > 1)
then 1161 self%mem(lb(1):ub(1),lb(2):ub(2),lb(3):lo(3),iv,it0,1) = exp( &
1162 log(self%mem(lb(1):ub(1),lb(2):ub(2),lb(3):lo(3),iv,it1,1))*2 - &
1163 log(self%mem(lb(1):ub(1),lb(2):ub(2),lb(3):lo(3),iv,it2,1)) )
1164 self%mem(lb(1):ub(1),lb(2):ub(2),uo(3):ub(3),iv,it0,1) = exp( &
1165 log(self%mem(lb(1):ub(1),lb(2):ub(2),uo(3):ub(3),iv,it1,1))*2 - &
1166 log(self%mem(lb(1):ub(1),lb(2):ub(2),uo(3):ub(3),iv,it2,1)) )
1171 if (self%gn(2) > 1)
then 1172 self%mem(lb(1):ub(1),lb(2):lo(2),li(3):ui(3),iv,it0,1) = exp( &
1173 log(self%mem(lb(1):ub(1),lb(2):lo(2),li(3):ui(3),iv,it1,1))*2 - &
1174 log(self%mem(lb(1):ub(1),lb(2):lo(2),li(3):ui(3),iv,it2,1)) )
1175 self%mem(lb(1):ub(1),uo(2):ub(2),li(3):ui(3),iv,it0,1) = exp( &
1176 log(self%mem(lb(1):ub(1),uo(2):ub(2),li(3):ui(3),iv,it1,1))*2 - &
1177 log(self%mem(lb(1):ub(1),uo(2):ub(2),li(3):ui(3),iv,it2,1)) )
1182 if (self%gn(1) > 1)
then 1183 self%mem(lb(1):lo(1),li(2):ui(2),li(3):ui(3),iv,it0,1) = exp( &
1184 log(self%mem(lb(1):lo(1),li(2):ui(2),li(3):ui(3),iv,it1,1))*2 - &
1185 log(self%mem(lb(1):lo(1),li(2):ui(2),li(3):ui(3),iv,it2,1)) )
1186 self%mem(uo(1):ub(1),li(2):ui(2),li(3):ui(3),iv,it0,1) = exp( &
1187 log(self%mem(uo(1):ub(1),li(2):ui(2),li(3):ui(3),iv,it1,1))*2 - &
1188 log(self%mem(uo(1):ub(1),li(2):ui(2),li(3):ui(3),iv,it2,1)) )
1194 if (self%gn(3) > 1)
then 1195 self%mem(lb(1):ub(1),lb(2):ub(2),lb(3):lo(3),iv,it0,1) = &
1196 self%mem(lb(1):ub(1),lb(2):ub(2),lb(3):lo(3),iv,it1,1)*2 - &
1197 self%mem(lb(1):ub(1),lb(2):ub(2),lb(3):lo(3),iv,it2,1)
1198 self%mem(lb(1):ub(1),lb(2):ub(2),uo(3):ub(3),iv,it0,1) = &
1199 self%mem(lb(1):ub(1),lb(2):ub(2),uo(3):ub(3),iv,it1,1)*2 - &
1200 self%mem(lb(1):ub(1),lb(2):ub(2),uo(3):ub(3),iv,it2,1)
1205 if (self%gn(2) > 1)
then 1206 self%mem(lb(1):ub(1),lb(2):lo(2),li(3):ui(3),iv,it0,1) = &
1207 self%mem(lb(1):ub(1),lb(2):lo(2),li(3):ui(3),iv,it1,1)*2 - &
1208 self%mem(lb(1):ub(1),lb(2):lo(2),li(3):ui(3),iv,it2,1)
1209 self%mem(lb(1):ub(1),uo(2):ub(2),li(3):ui(3),iv,it0,1) = &
1210 self%mem(lb(1):ub(1),uo(2):ub(2),li(3):ui(3),iv,it1,1)*2 - &
1211 self%mem(lb(1):ub(1),uo(2):ub(2),li(3):ui(3),iv,it2,1)
1216 if (self%gn(1) > 1)
then 1217 self%mem(lb(1):lo(1),li(2):ui(2),li(3):ui(3),iv,it0,1) = &
1218 self%mem(lb(1):lo(1),li(2):ui(2),li(3):ui(3),iv,it1,1)*2 - &
1219 self%mem(lb(1):lo(1),li(2):ui(2),li(3):ui(3),iv,it2,1)
1220 self%mem(uo(1):ub(1),li(2):ui(2),li(3):ui(3),iv,it0,1) = &
1221 self%mem(uo(1):ub(1),li(2):ui(2),li(3):ui(3),iv,it1,1)*2 - &
1222 self%mem(uo(1):ub(1),li(2):ui(2),li(3):ui(3),iv,it2,1)
1226 call trace%end (itimer)
1227 END SUBROUTINE extrapolate_guards
1232 SUBROUTINE allocate_mesg (self)
1235 allocate (self%mesg)
1236 self%mesg%nbuf = product(self%gn)
1237 if (kind(self%mem)==8) &
1238 self%mesg%nbuf = self%mesg%nbuf * 2
1239 self%mesg%nbuf = n_header + self%nv*self%mesg%nbuf
1240 allocate (self%mesg%buffer(self%mesg%nbuf))
1241 call io%bits_mem (storage_size(self%mesg%buffer), product(shape(self%mesg%buffer)))
1242 END SUBROUTINE allocate_mesg
1247 SUBROUTINE pack (self, mesg)
1249 class(
mesg_t),
pointer:: mesg
1252 integer:: n_buf, ibuf, iv, it
1254 integer,
save:: itimer=0
1258 call trace%begin (
'patch_t%pack', itimer=itimer)
1260 n_buf = product(self%gn)
1261 if (kind(self%mem)==8) &
1263 mesg%nbuf = n_header + self%nv*n_buf
1264 if (self%is_set(bits%swap_request) .and. self%is_set(bits%boundary))
then 1265 mesg%nbuf = mesg%nbuf+ self%nv*n_buf
1267 allocate (mesg%buffer(mesg%nbuf))
1268 if (mpi_mesg%nbuf==0)
then 1270 if (mpi_mesg%nbuf==0)
then 1271 mpi_mesg%nbuf = mesg%nbuf
1272 write (io_unit%log,*) &
1273 'patch_t%pack: setting mpi_mesg%nbuf =', mpi_mesg%nbuf
1277 call io%bits_mem (storage_size(mesg%buffer), product(shape(mesg%buffer)))
1278 allocate (mesg%reqs(self%n_nbors))
1284 if (self%debug(2))
then 1285 write (io_unit%log,1) wallclock(), &
1286 'DBG patch_t%pack, id, nq, time, status:', &
1287 self%id, self%nq, self%time, &
1288 self%is_set (bits%internal), &
1289 self%is_set (bits%boundary), &
1290 self%is_set (bits%virtual), &
1291 self%is_set (bits%external), &
1292 self%is_set (bits%swap_request)
1293 1
format(f12.6,2x,a,2i9,f12.6,2x,5l1)
1296 call patch_to_header (self, header)
1298 call anonymous_copy (n_header, header, mesg%buffer(ibuf))
1299 ibuf = ibuf + n_header
1305 call anonymous_copy (n_buf, self%mem(:,:,:,iv,self%it,1), mesg%buffer(ibuf))
1308 if (self%is_set(bits%swap_request) .and. self%is_set(bits%boundary))
then 1309 it = mod(self%it-1+self%nt-1,self%nt)+1
1310 if (io%verbose>1)
then 1311 write (io_unit%log,*) self%id,
'extra time slot packed', it
1315 call anonymous_copy (n_buf, self%mem(:,:,:,iv,it,1), mesg%buffer(ibuf))
1319 if (io%verbose > 2) &
1320 write (io_unit%log,*) self%id,
'after pack', minval(self%mem(:,:,:,1,self%it,1))
1321 call trace%end (itimer)
1327 INTEGER FUNCTION unpack_id (mesg)
1328 class(
mesg_t),
pointer:: mesg
1332 call trace%begin (
'patch_t%unpack_id')
1333 call anonymous_copy (n_header, mesg%buffer, header)
1334 unpack_id = header%id
1336 END FUNCTION unpack_id
1341 SUBROUTINE unpack (self, mesg)
1343 class(
mesg_t),
pointer:: mesg
1346 integer:: n_buf, ibuf, iv, it
1347 integer,
dimension(:),
pointer:: buffer
1348 integer,
save:: itimer=0
1352 call trace%begin (
'patch_t%unpack', itimer=itimer)
1353 call self%set (bits%busy)
1354 self%unpack_time = wallclock()
1355 buffer => mesg%buffer
1356 n_buf = product(self%gn)
1357 if (kind(self%mem)==8) &
1363 call anonymous_copy (n_header, buffer(ibuf), header)
1364 ibuf = ibuf + n_header
1368 call self%header_to_patch (header)
1369 if (io%verbose > 1) &
1370 write (io_unit%log,
'(f12.6,3x,a,i7,2i4,2x,3l1,2i4,f8.4,f11.6)') &
1371 wallclock(),
'header_to_patch: id, nq, BVS, nsnd, nrcv, latency, time =', &
1372 self%id, mesg%sender, self%nq, &
1373 self%is_set(bits%boundary), &
1374 self%is_set(bits%virtual), &
1375 self%is_set(bits%swap_request), &
1376 mpi_mesg%sent_list%n, mpi_mesg%recv_list%n, self%latency, self%time
1377 self%wc_last = wallclock()
1378 if (io%verbose > 1)
write (io_unit%log,1)
'patch_mod::unpk id, nq, time after:', &
1379 self%id, self%nq, self%time, &
1380 self%is_set (bits%internal), &
1381 self%is_set (bits%boundary), &
1382 self%is_set (bits%virtual), &
1383 self%is_set (bits%external), &
1384 self%is_set (bits%swap_request), &
1386 1
format(a,2i9,f12.6,2x,5l1,f12.1)
1391 call self%mem_lock(self%it)%set (
'unpack')
1393 call anonymous_copy (n_buf, buffer(ibuf), self%mem(:,:,:,iv,self%it,1))
1399 call self%mem_lock(self%it)%unset (
'unpack')
1400 if (self%is_set(bits%swap_request) .and. self%is_set(bits%virtual))
then 1401 it = mod(self%it-1+self%nt-1,self%nt)+1
1402 call self%mem_lock(it)%set (
'unpack')
1403 if (io%verbose>0)
then 1404 write (io_unit%log,*) self%id,
'extra time slot unpacked', it
1408 call anonymous_copy (n_buf, buffer(ibuf), self%mem(:,:,:,iv,it,1))
1411 call self%mem_lock(it)%unset (
'unpack')
1413 if (io%verbose > 2) &
1414 write (io_unit%log,*) self%id,
'after unpack', minval(self%mem(:,:,:,1,self%it,1))
1415 call self%clear (bits%busy)
1416 call trace%end (itimer)
1417 END SUBROUTINE unpack
1426 SUBROUTINE info (self, nq, ntask, experiment_name)
1428 integer,
optional:: nq, ntask
1429 character(len=64),
optional:: experiment_name
1431 integer:: iv, nv, iprint_l
1432 integer:: it0, it1, it
1434 real,
dimension(:,:,:),
pointer:: df
1435 integer,
parameter:: max_lines=50
1436 integer,
save:: counts(6)=0, id_prv=0
1437 integer:: counts_l(6), i
1438 character(len=120):: fmt
1439 real(8):: time, print_next
1441 type(
lock_t),
save:: lock
1442 integer,
save:: itimer=0
1444 call trace%begin (
'patch_t%info', itimer=itimer)
1445 if (io_unit%do_validate)
then 1446 it0 = merge(self%iit(self%nt-2), 1, self%nt>2)
1447 it1 = merge(self%iit(self%nt-1), 1, self%nt>1)
1448 if (self%nw==1)
then 1449 allocate (df(
size(self%mem,1),
size(self%mem,2),
size(self%mem,3)))
1452 if (self%nw==1)
then 1453 df = self%mem(:,:,:,iv,it1,1) &
1454 - self%mem(:,:,:,iv,it0,1)
1455 fmin = self%fminval(df)
1456 fmax = self%fmaxval(df)
1458 fmin = self%fminval(self%mem(:,:,:,iv,it0,self%nw))
1459 fmax = self%fmaxval(self%mem(:,:,:,iv,it0,self%nw))
1461 write(io%output,2) &
1462 'upd',self%id, self%level, iv, self%istep, self%time-self%dtime, &
1463 self%dtime, self%u_max, fmin, fmax, &
1464 self%fminval(self%mem(:,:,:,iv,it0,1)), &
1465 self%fmaxval(self%mem(:,:,:,iv,it0,1)), &
1466 real(timer%n_update), ntask, nq, self%is_set(bits%boundary), self%is_set(bits%virtual)
1467 2
format(a,i8,2i4,i3,1p,e14.5,e12.4,0p,f7.2,2x,1p,2e12.4,1x,2e12.4,2x,e10.2,2i5,l3,l1,o12)
1469 if (self%nw==1)
then 1472 call trace%end (itimer)
1475 1
format(a,i8,2i4,i3,f12.6,1p,e12.4,0p,f7.2,2x,1p,2e10.2,2x,2e11.3,2x,e10.2,9i5,2i4,l3,l1,o12)
1476 if (self%time < 1d4)
then 1477 fmt =
'(a,i8,2i4,i3,f12.6,1p,e12.4,0p,f7.2,2x,1p,2e10.2,2x,2e11.3,2x,e10.2,i6,3i5,i6,4i5,2i4,l3,l1,o12)' 1479 fmt =
'(a,i8,2i4,i3,1p,2e12.5, 0p,f7.2,2x,1p,2e10.2,2x,2e11.3,2x,e10.2,i6,3i5,i6,4i5,2i4,l3,l1,o12)' 1482 if ((io%verbose>=0).and.(self%is_clear(bits%no_density)))
then 1484 io%dmin = min(io%dmin,
real(self%fminval(self%mem(:,:,:,self%idx%d,self%it,1)),kind=4))
1486 io%dmax = max(io%dmax,
real(self%fmaxval(self%mem(:,:,:,self%idx%d,self%it,1)),kind=4))
1490 if (self%id == id_prv)
then 1497 if (print_every > 0)
then 1500 if (iprint_l >= print_every .and. print_every > 0)
then 1513 else if (io%print_time > 0.0)
then 1514 if (print_every == 0 .and. self%time >= io%print_next)
then 1515 call lock%set (
'print')
1516 if (self%time >= io%print_next)
then 1518 print_next = (nint(self%time/io%print_time)+1)*io%print_time
1519 io%print_next = print_next
1521 call lock%unset (
'print')
1526 else if (self%id==io%id)
then 1529 if (print_it .or. io%verbose>0)
then 1531 if (io%verbose>2) nv = self%nv
1534 [mpi_mesg%n_check, mpi_mesg%n_ready, mpi_mesg%n_update, &
1535 mpi_mesg%n_send , mpi_mesg%n_recv , mpi_mesg%n_unpk ]
1538 counts(i) = counts(i) + counts_l(i)
1541 timer%n_lines = timer%n_lines-1
1542 if (timer%n_lines==0)
then 1545 timer%n_lines = max_lines
1547 write(io%output,
'(a)') &
1548 '0....+....1....+....2....+....3....+....4....+....5....+....' &
1549 //
'6....+....7....+....8....+....9....+....0....+....1 ntsk nq ' &
1550 //
' chk rdy upd snt prb unp snq rcq unq BV' 1559 if (io%print_time == 0.0 .and. io%verbose > 0)
then 1560 time = self%time - self%dtime
1565 write(io%output,fmt) &
1566 'upd',self%id, self%level, iv, omp_mythread, time, &
1567 self%dtime, log10(max(self%u_max,1e-30)), &
1569 self%fminval(self%mem(:,:,:,iv,self%it ,1)), &
1570 self%fmaxval(self%mem(:,:,:,iv,self%it ,1)), &
1571 real(timer%n_update), ntask, nq, counts, &
1572 mpi_mesg%nq_send, mpi_mesg%nq_recv, mpi_mesg%unpk_list%n, &
1573 self%is_set(bits%boundary), self%is_set(bits%virtual)
1581 io%dtime = huge(1d0)
1588 mpi_mesg%n_check=0; mpi_mesg%n_ready=0; mpi_mesg%n_update=0
1589 mpi_mesg%n_send=0 ; mpi_mesg%n_recv=0 ; mpi_mesg%n_unpk=0
1598 if (self%dtime < io%dtime .and. id_prv /= 0)
then 1602 io%dtime = self%dtime
1605 call trace%end (itimer)
1611 LOGICAL FUNCTION contains1 (self, point)
1614 contains1 = all(abs(point-self%position) < 0.5_8*self%size)
1615 END FUNCTION contains1
1620 LOGICAL FUNCTION contains2 (self, patch)
1622 contains2 = contains1(self, patch%position)
1623 END FUNCTION contains2
1631 FUNCTION overlaps (self, task)
1635 real(8):: distance, box, limit
1643 distance = self%centre_nat(i) - task%centre_nat(i)
1645 if (self%periodic(i))
then 1646 distance = modulo(distance+0.5_8*box, box) - 0.5_8*box
1648 limit = 0.5_8*(self%size(i)+task%size(i)+self%ds(i)+task%ds(i))
1649 if (abs(distance) > limit)
then 1655 overlaps = self%task_t%overlaps(task)
1658 END FUNCTION overlaps
1666 FUNCTION distance_curvilinear (self, task)
RESULT (out)
1671 call trace%begin (
'patch_t%distance_curvilinear',2)
1675 out = self%centre_nat-task%centre_nat
1678 out = self%position-task%position
1682 where (self%periodic .and. self%box > 0d0)
1683 out = modulo(out+0.5_8*self%box, self%box) - 0.5_8*self%box
1685 if (io%verbose>1 .and. (task%id==io%id_debug .or. self%id==io%id_debug))
then 1686 write(io%output,1) &
1687 'distance: self', self%id, self%position, self%size, self%box
1688 write(io%output,1) &
1689 'distance: task', task%id, task%position, task%size, task%box, out
1690 1
format(a,i5,4(3x,3f10.2))
1693 END FUNCTION distance_curvilinear
1703 FUNCTION index (self, other, oindex, f)
1705 real,
dimension(3):: position, f, pos, pos_lb
1706 integer,
dimension(3):: index, oindex
1708 pos = other%position + (oindex-other%offset)*other%ds
1709 pos = pos-self%position
1710 pos = modulo(pos+0.5*self%box,self%box) - 0.5*self%box
1711 f = self%offset + pos/self%ds
1720 FUNCTION nearest (self, other, index)
RESULT (out)
1722 integer,
dimension(3):: out, index
1723 real,
dimension(3):: position, pos_lb, f, pos
1725 pos = other%position + (index-other%offset)*other%ds
1726 pos = pos-self%position
1727 pos = modulo(pos+0.5*self%box,self%box) - 0.5*self%box
1728 f = self%offset + pos/self%ds
1730 END FUNCTION nearest
1735 FUNCTION lowest (self, other, index)
RESULT (out)
1737 integer,
dimension(3):: out, index
1738 real,
dimension(3):: position, pos_lb, f, pos
1740 pos = other%position + (index-other%offset-0.49)*other%ds
1741 pos = pos-self%position
1742 pos = modulo(pos+0.5*self%box,self%box) - 0.5*self%box
1743 f = self%offset + pos/self%ds
1749 FUNCTION highest (self, other, index)
RESULT (out)
1751 integer,
dimension(3):: out, index
1752 real,
dimension(3):: position, pos_lb, f, pos
1754 pos = other%position + (index-other%offset+0.49)*other%ds
1755 pos = pos-self%position
1756 pos = modulo(pos+0.5*self%box,self%box) - 0.5*self%box
1757 f = self%offset + pos/self%ds
1759 END FUNCTION highest
1764 SUBROUTINE courant_condition (self, detailed_timer)
1766 logical,
optional:: detailed_timer
1767 real(8):: dtime, ds(3), courant
1769 integer,
save:: itimer=0
1770 logical,
parameter:: save_dbg=.false.
1772 call trace%begin(
'patch_t%courant_condition', itimer=itimer)
1774 write (io_unit%log,*) self%id,
' courant: u_max, ds', self%u_max, minval(self%ds)
1776 ds(2) = self%ds(2)*self%mesh(1)%h2c(self%mesh(1)%li)
1777 ds(3) = self%ds(3)*self%mesh(1)%h31c(self%mesh(1)%li)*abs(self%mesh(2)%h32c(self%mesh(2)%li))
1778 if (trim(self%mesh(3)%type) ==
'spherical_mesh')
then 1779 do i2=self%mesh(2)%li,self%mesh(2)%ui
1780 ds(3) = min(ds(3),self%ds(3)*self%mesh(1)%h31c(i2)*abs(self%mesh(2)%h32c(i2)))
1783 if (self%dt_fixed > 0.0)
then 1784 dtime = self%dt_fixed
1785 courant = dtime*self%u_max/minval(ds)
1786 else if (self%u_fixed > 0.0)
then 1787 dtime = self%courant*minval(ds)/self%u_fixed
1788 courant = self%courant
1790 dtime = self%courant*minval(ds)/self%u_max*self%safety_factor
1791 courant = self%courant
1795 if (self%nt > 1 .and. self%istep > 5)
then 1796 if (dtime < 0.5*self%dt(self%iit(self%nt-1)) .and. self%dt_fixed < 0.0)
then 1797 self%safety_factor = 0.5
1799 if (self%safety_factor < 0.6 .and. save_dbg)
then 1801 if (self%safety_factor == 0.5)
then 1802 write (io_unit%dbg) self%id, self%gn, self%ipos
1803 write (io_unit%dbg) self%t(self%nt-3), self%mem(:,:,:,:,self%iit(self%nt-3),1)
1804 write (io_unit%dbg) self%id, self%gn, self%ipos
1805 write (io_unit%dbg) self%t(self%nt-2), self%mem(:,:,:,:,self%iit(self%nt-2),1)
1807 write (io_unit%dbg) self%id, self%gn, self%ipos
1808 write (io_unit%dbg) self%time, self%mem(:,:,:,:,self%it,1)
1810 if (io%verbose>2 .or. (self%istep>5)) &
1811 write(io%output,
'(a,i6,i3,f12.6,1p,3g11.2,0p,1x,3f12.6,1x,1p,7g12.3)') &
1812 ' courant: ', self%id, self%it, &
1813 self%time, dtime, courant, self%u_max, self%position, self%dt(self%iit)
1816 self%safety_factor = self%safety_factor**0.9
1817 self%get_umax_location = (self%safety_factor < 0.9)
1822 if (self%time < self%sync_time .and. self%time+self%dtime > self%sync_time)
then 1823 self%dtime = self%sync_time - self%time
1824 self%syncing = .true.
1826 write (io_unit%mpi,*)
'patch_t%courant_condition: syncing turned on for', &
1832 if (io%do_output .and. io%out_time > 0.0) &
1833 self%dtime = min(self%dtime, io%out_time)
1834 self%dt(self%it) = self%dtime
1836 write(io%output,
'(a,i7,f12.6,f8.4,1p,e10.2)')
'courant:', &
1837 self%id, self%dtime, self%courant, self%u_max
1838 call trace%end (itimer)
1839 END SUBROUTINE courant_condition
1842 FUNCTION fsum4 (self, f)
RESULT (out)
1845 real,
dimension(:,:,:),
intent(in):: f
1848 if (self%is_periodic())
then 1850 out = sum(
real(f(self%mesh(1)%li:self%mesh(1)%li+n(1), &
self%mesh(2)%li:self%mesh(2)%li+n(2), &
self%mesh(3)%li:self%mesh(3)%li+n(3)),kind=8))
1852 out = sum(
real(f(self%mesh(1)%li:self%mesh(1)%ui, &
self%mesh(2)%li:self%mesh(2)%ui, &
self%mesh(3)%li:self%mesh(3)%ui),kind=8))
1856 FUNCTION fsum8 (self, f)
RESULT (out)
1859 real(8),
dimension(:,:,:),
intent(in):: f
1862 if (self%is_periodic())
then 1864 out = sum(
real(f(self%mesh(1)%li:self%mesh(1)%li+n(1), &
self%mesh(2)%li:self%mesh(2)%li+n(2), &
self%mesh(3)%li:self%mesh(3)%li+n(3)),kind=8))
1866 out = sum(
real(f(self%mesh(1)%li:self%mesh(1)%ui, &
self%mesh(2)%li:self%mesh(2)%ui, &
self%mesh(3)%li:self%mesh(3)%ui),kind=8))
1871 FUNCTION faver4 (self, f)
RESULT (out)
1874 real,
dimension(:,:,:),
intent(in):: f
1876 out = self%fsum4(f)/product(self%ncell)
1879 FUNCTION faver8 (self, f)
RESULT (out)
1882 real(8),
dimension(:,:,:),
intent(in):: f
1884 out = self%fsum8(f)/product(self%ncell)
1888 FUNCTION fminvali (self, iv)
RESULT (fminval)
1893 fminval = self%fminvalvar (self%mem(:,:,:,iv,self%it,1))
1894 END FUNCTION fminvali
1896 FUNCTION fmaxvali (self, iv)
RESULT (fmaxval)
1901 fmaxval = self%fmaxvalvar (self%mem(:,:,:,iv,self%it,1))
1902 END FUNCTION fmaxvali
1904 FUNCTION fmaxval4 (self, f, outer)
RESULT (fmaxval)
1907 real,
dimension(:,:,:),
intent(in):: f
1908 logical,
optional:: outer
1911 if (
present(outer))
then 1916 if (self%boundaries%is_set (bits%spherical))
then 1917 fmaxval = maxval(f, mask=self%boundaries%mask(self%mesh))
1918 else if (outer_l)
then 1919 fmaxval = maxval(f(self%mesh(1)%lo:self%mesh(1)%uo, &
1920 self%mesh(2)%lo:self%mesh(2)%uo, &
1921 self%mesh(3)%lo:self%mesh(3)%uo))
1923 fmaxval = maxval(f(self%mesh(1)%li:self%mesh(1)%ui, &
1924 self%mesh(2)%li:self%mesh(2)%ui, &
1925 self%mesh(3)%li:self%mesh(3)%ui))
1927 END FUNCTION fmaxval4
1929 FUNCTION fminval4 (self, f)
RESULT(fminval)
1932 real,
dimension(:,:,:),
intent(in):: f
1934 if (self%boundaries%is_set (bits%spherical))
then 1935 fminval = minval(f, mask=self%boundaries%mask(self%mesh))
1937 fminval = minval(f(self%mesh(1)%li:self%mesh(1)%ui, &
1938 self%mesh(2)%li:self%mesh(2)%ui, &
1939 self%mesh(3)%li:self%mesh(3)%ui))
1941 END FUNCTION fminval4
1943 FUNCTION fmaxval8 (self, f, outer)
RESULT (fmaxval)
1946 real(8),
dimension(:,:,:),
intent(in):: f
1947 logical,
optional:: outer
1950 if (
present(outer))
then 1955 if (self%boundaries%is_set (bits%spherical))
then 1956 fmaxval = maxval(f, mask=self%boundaries%mask(self%mesh))
1957 else if (outer_l)
then 1958 fmaxval = maxval(f(self%mesh(1)%lo:self%mesh(1)%uo, &
1959 self%mesh(2)%lo:self%mesh(2)%uo, &
1960 self%mesh(3)%lo:self%mesh(3)%uo))
1962 fmaxval = maxval(f(self%mesh(1)%li:self%mesh(1)%ui, &
1963 self%mesh(2)%li:self%mesh(2)%ui, &
1964 self%mesh(3)%li:self%mesh(3)%ui))
1966 END FUNCTION fmaxval8
1968 FUNCTION fminval8 (self, f)
RESULT(fminval)
1971 real(8),
dimension(:,:,:),
intent(in):: f
1973 if (self%boundaries%is_set (bits%spherical))
then 1974 fminval = minval(f, mask=self%boundaries%mask(self%mesh))
1976 fminval = minval(f(self%mesh(1)%li:self%mesh(1)%ui, &
1977 self%mesh(2)%li:self%mesh(2)%ui, &
1978 self%mesh(3)%li:self%mesh(3)%ui))
1980 END FUNCTION fminval8
1985 SUBROUTINE check_nan (self, panic, label, always)
1988 character(len=*),
optional:: label
1989 logical,
optional:: always
1991 integer:: ix, iy, iz, iv, ii(3)
1993 integer,
save:: itimer=0
1995 if (.not. do_check_nan)
then 1996 if (io%verbose < 2 .and. .not.
present(always))
return 1998 if (.not.
associated(self%mem))
return 1999 call timer%begin (
'patch_t%check_nan', itimer=itimer)
2000 if (io%verbose>1)
then 2001 if (
present(label))
then 2002 write(io%output,*) self%id, trim(label)//
', min(d), max(d):', &
2003 minval(self%mem(:,:,:,self%idx%d,self%it,1)), &
2004 maxval(self%mem(:,:,:,self%idx%d,self%it,1))
2006 write(io%output,*) self%id,
'min(d), max(d):', &
2007 minval(self%mem(:,:,:,self%idx%d,self%it,1)), &
2008 maxval(self%mem(:,:,:,self%idx%d,self%it,1))
2016 if (isinf(self%mem(ix,iy,iz,iv,self%it,1)))
then 2017 if (io%verbose>0)
then 2018 if (
present(label))
then 2019 write(io%output,1) trim(label)//
' ERROR: NaN id,ix,iy,iz,iv:', &
2021 1
format(a,i6,3i4,i5)
2023 write(io%output,1)
'ERROR: NaN ix,iy,iz,iv:', ix,iy,iz,iv
2037 if (
present(label))
then 2038 write(io%output,*) self%id, trim(label)//
' ERROR NaN, min(d), max(d):', &
2039 minval(self%mem(:,:,:,self%idx%d,self%it,1)), &
2040 maxval(self%mem(:,:,:,self%idx%d,self%it,1)), ii
2042 write(io%output,*) self%id,
' ERROR NaN, min(d), max(d):', &
2043 minval(self%mem(:,:,:,self%idx%d,self%it,1)), &
2044 maxval(self%mem(:,:,:,self%idx%d,self%it,1)), ii
2046 write(io%output,*) self%unsigned
2048 call io%abort (
'found NaN -- check rank & thread logs!')
2050 call timer%end (itimer)
2051 END SUBROUTINE check_nan
2053 LOGICAL FUNCTION isinf(a)
2054 use ieee_arithmetic
, only : ieee_is_nan
2055 real(kind=KindScalarVar):: a
2056 isinf = ieee_is_nan(a) .or. (abs(a) > huge(1.0_kindscalarvar))
2063 SUBROUTINE check_density (self, label)
2065 character(len=*),
optional:: label
2066 real(kind=KindScalarVar):: mind, posd
2067 real:: tiny = 1.0e-32
2068 integer:: ix, iy, iz, iv, ii(3), ntot, nins, l(3), u(3), it, jt, iw
2069 logical:: inside, nan, panic
2070 integer,
save:: itimer=0
2072 if ((io%verbose < 2 .or. io%print_time > 0.0 .or. print_every > 1) .and. &
2073 .not.do_check_nan)
return 2074 if (self%is_set(bits%no_density+bits%frozen))
return 2075 if (.not.
associated(self%mem))
return 2076 call timer%begin (
'patch_t%check_density', itimer=itimer)
2079 mind = minval(self%mem(l(1):u(1),l(2):u(2),l(3):u(3),self%idx%d,self%it,1))
2080 if (
present(label) .and. io%verbose>1)
write (io_unit%log,*) label, self%id, mind
2082 if (mind <= 0.0)
then 2095 inside = all(ii>self%mesh%lo .and. ii<self%mesh%uo)
2096 posd = merge(posd, max(posd,self%mem(ix,iy,iz,self%idx%d,self%it,1)), inside)
2097 ntot = ntot + merge(1,0,self%mem(ix,iy,iz,self%idx%d,self%it,1)<=0.0)
2098 nins = nins + merge(1,0,self%mem(ix,iy,iz,self%idx%d,self%it,1)<=0.0 .and. inside)
2108 1
format(a,i8,3x,a,i8,a,i8,a)
2109 if (
present(label))
then 2111 write (io_unit%output,1)
'ID ',self%id, &
2112 'ERROR: density non-positive in',ntot,
' points,',nins,
' inside, '// &
2114 else if (io%verbose>0)
then 2115 write (io_unit%output,1)
'ID ',self%id, &
2116 'WARNING: density non-positive in',ntot,
' points,',nins,
' inside, '// &
2121 write (io_unit%output,1)
'ID ',self%id, &
2122 'WARNING: density non-positive in',ntot,
' points,',nins,
' inside' 2123 else if (io%verbose>0)
then 2124 write (io_unit%output,1)
'ID ',self%id, &
2125 'WARNING: density non-positive in',ntot,
' points,',nins,
' inside' 2131 call self%check_nan (panic, label)
2137 write (io_unit%dump) self%id, self%gn, self%nv, self%nt, self%nw, self%istep
2138 write (io_unit%dump) self%time, self%dtime
2144 write (io_unit%dump) self%mem(:,:,:,:,jt,iw)
2148 if (nins > 0 .or. io%verbose > 0) &
2149 print *, mpi%rank, self%id, &
2150 'WARNING: found non-positive density points -- check rank and thread logs' 2152 write (stderr,*)
'check dump, rank and thread logs for rank', mpi%rank
2156 flush (io_unit%dump)
2157 call io%abort (
'found non-positive internal density points')
2160 call timer%end (itimer)
2161 END SUBROUTINE check_density
2166 SUBROUTINE stats_patch (self, label)
2168 character(len=*):: label
2170 real(8),
dimension(4):: si, sb
2172 call scalar_stats (self%mem(:,:,:,iv,self%it,1), self%mesh%li, self%mesh%ui, si)
2173 call scalar_stats (self%mem(:,:,:,iv,self%it,1), self%mesh%lb, self%mesh%ub, sb)
2174 write (io_unit%log,
'(i7,2x,a,i3,1p,2(3x,4g12.4))') &
2175 self%id, trim(label)//
' stats: iv, min, aver, rms, max =', iv, si, sb
2177 END SUBROUTINE stats_patch
2182 SUBROUTINE stats_scalar (self, var, label)
2184 real(kind=KindScalarVar),
dimension(:,:,:):: var
2185 character(len=*):: label
2187 real(8),
dimension(4):: si, sb
2189 call scalar_stats (var, self%mesh%li, self%mesh%ui, si)
2190 call scalar_stats (var, self%mesh%lb, self%mesh%ub, sb)
2191 write (io_unit%output,
'(i7,2x,a,1p,2(3x,4g12.4))') &
2192 self%id, trim(label)//
' stats: min, aver, rms, max =', si, sb
2193 END SUBROUTINE stats_scalar
2198 SUBROUTINE stats_vector (self, var, label)
2200 real(kind=KindScalarVar),
dimension(:,:,:,:):: var
2201 character(len=*):: label
2204 real(8),
dimension(4):: si, sb
2207 call scalar_stats (var(:,:,:,iv), self%mesh%li, self%mesh%ui, si)
2208 call scalar_stats (var(:,:,:,iv), self%mesh%lb, self%mesh%ub, sb)
2209 write (io_unit%output,
'(i7,2x,a,i3,1p,2(3x,4g12.4))') &
2210 self%id, trim(label)//
' stats: iv, min, aver, rms, max =', iv, si, sb
2212 END SUBROUTINE stats_vector
2216 SUBROUTINE scalar_stats (f, l, u, s)
2217 real(kind=KindScalarVar),
dimension(:,:,:):: f
2218 integer:: l(3), u(3)
2219 real(8),
dimension(4):: s
2220 integer:: ix, iy, iz
2228 s(2) = s(2) + f(ix,iy,iz)
2229 s(3) = s(3) + f(ix,iy,iz)**2
2230 s(1) = min(s(1),f(ix,iy,iz))
2231 s(4) = max(s(4),f(ix,iy,iz))
2235 s(2) = s(2)/product(u-l+1)
2236 s(3) = sqrt(max(0d0,s(3)/product(u-l+1)-s(2)**2))
2237 END SUBROUTINE scalar_stats
2243 SUBROUTINE mapvectors (self, in, target, ii, jj, out)
2245 class(
patch_t),
pointer :: target
2246 real(kind=KindScalarVar),
dimension(:),
intent(in):: in
2247 real(kind=KindScalarVar),
intent(out):: out(size(in))
2248 integer,
intent(in) :: ii(3), jj(3)
2250 call trace%begin(
'patch_t%MapVectors', 2)
2252 select case (self%mesh_type)
2254 if (
target%mesh_type == mesh_types%Cartesian)
then 2258 call mpi%abort(
'Not implemented error! patch_t%MapVectors')
2261 call mpi%abort(
'Not implemented error! patch_t%MapVectors')
2265 END SUBROUTINE mapvectors
2268 FUNCTION beyond_patch_edge (self, idir, position)
result (outside)
2270 integer,
intent(in) :: idir
2271 real(8),
intent(in) :: position(3)
2273 real(8) :: edge, curvpos(3)
2276 if (self%n(idir) == 1)
then 2281 select case (self%mesh_type)
2283 edge = 0.5 * self%gsize(idir)
2284 outside = abs(position(idir)-self%position(idir)) > edge
2286 curvpos = position - self%centre_cart
2287 curvpos = cartesiantospherical(curvpos(1), curvpos(2), curvpos(3))
2290 edge = self%gsize(idir) + self%llc_nat(idir)
2291 outside = abs(curvpos(idir)) > edge .or. abs(curvpos(idir)) < self%llc_nat(idir)
2293 if (
associated(self%mesh))
then 2294 edge = self%gsize(idir) + self%mesh(2)%xi
2296 edge = self%gsize(idir)
2298 outside = abs(curvpos(idir)) > edge
2300 if (
associated(self%mesh))
then 2301 edge = self%gsize(idir) + self%mesh(3)%zeta
2303 edge = self%gsize(idir)
2305 outside = abs(curvpos(idir)) > edge
2307 outside = abs(curvpos(idir)) > edge
2309 curvpos = position - self%centre_cart
2310 curvpos = cartesiantocylindrical(curvpos(1), curvpos(2), curvpos(3))
2313 edge = 0.5 * self%gsize(idir)
2314 outside = abs(curvpos(idir)) > edge
2316 edge = self%gsize(idir) + self%llc_nat(idir)
2317 outside = abs(curvpos(idir)) > edge .or. abs(curvpos(idir)) < self%llc_nat(idir)
2319 if (
associated(self%mesh))
then 2320 edge = self%gsize(idir) + self%mesh(3)%zeta
2322 edge = self%gsize(idir)
2324 outside = abs(curvpos(idir)) > edge
2327 END FUNCTION beyond_patch_edge
2333 SUBROUTINE log_interpolate (self, time, i, a)
2337 real(kind=KindScalarVar),
dimension(:,:,:),
pointer :: a, b
2341 integer,
save :: itimer=0
2343 call trace%begin (
'patch_t%log_interpolate', itimer=itimer)
2344 call self%time_interval (time, jt, pt)
2345 if (io%verbose > 1)
then 2346 b => self%mem(:,:,:,i,jt(1),1)
2347 write(io%output,*) self%id,
'log_interp(1)', jt(1), minval(b), maxval(b)
2348 b => self%mem(:,:,:,i,jt(2),1)
2349 write(io%output,*) self%id,
'log_interp(2)', jt(2), minval(b), maxval(b)
2352 a = log(self%mem(:,:,:,i,jt(1),1))*pt(1) &
2353 + log(self%mem(:,:,:,i,jt(2),1))*pt(2)
2354 call trace%end (itimer)
2355 END SUBROUTINE log_interpolate
2361 SUBROUTINE time_interval (self, time, jt, pt, all)
2366 logical,
optional :: all
2368 integer :: i, j, n, iit(self%nt-1)
2369 real(8) :: t(self%nt-1), dt
2370 integer,
save :: itimer=0
2373 if (self%nt==1)
then 2377 call self%timeslots (iit, t)
2388 if (t(i) < time) j=i
2395 pt(2) = (time-t(j))/dt
2399 if ((io%verbose>3 .or. io_unit%do_validate) .and. pt(2)>0.0)
then 2400 write(io%output,
'(a,i6,i6,1p,e14.6,i3,0p,2f8.4,2x,1p,7e14.6)') &
2401 'time_interval: times =', self%id, self%istep, time, j, pt, t
2404 END SUBROUTINE time_interval
2408 FUNCTION index_stagger (self, iv)
RESULT (out)
2415 out(i) = self%mesh(i)%h(iv)
2417 END FUNCTION index_stagger
2428 SUBROUTINE index_space (self, pos, iv, i, p)
2436 real,
parameter :: eps=0.001
2437 integer,
save :: itimer=0
2440 h = self%index_stagger(iv)
2441 p = pos - self%position
2442 p = modulo(p+0.5_8*self%mesh%b,self%mesh%b) &
2444 p = p/self%mesh%d + self%mesh%o - h
2447 i = max(self%mesh%li, &
2448 min(self%mesh%ui-1,i))
2453 if (zero_order_extrap .and. .not.io_unit%do_validate) &
2454 p = min(p,max(p,-0.01),1.01)
2456 END SUBROUTINE index_space
2462 SUBROUTINE index_space_of (self, target, ii, iv, jj, pp)
2464 class(
patch_t),
pointer :: target
2470 class(
mesh_t),
pointer :: mt, ms
2473 real,
parameter :: eps=0.001
2474 integer,
save :: itimer=0
2477 mt =>
target%mesh(i)
2479 p = mt%p - ms%p + (ii(i) - mt%o + mt%h(iv))*mt%d
2480 p = p/ms%d + ms%o - ms%h(iv) + eps
2485 if (zero_order_extrap .and. .not.io_unit%do_validate) &
2486 p = min(max(p,-0.01),1.01)
2490 END SUBROUTINE index_space_of
2498 FUNCTION index_only (self, pos, iv, roundup, nowrap)
RESULT (ii)
2500 real(8),
intent(in) :: pos(3)
2501 integer,
optional :: iv
2502 logical,
optional :: roundup, nowrap
2506 real(8) :: p(3), h(3)
2507 class(
mesh_t),
pointer :: mesh
2510 h = self%index_stagger(iv)
2512 mesh => self%mesh(i)
2513 p(i) = pos(i) - self%position(i)
2514 if (self%periodic(i))
then 2515 if (.not.
present(nowrap)) &
2516 p(i) = modulo(p(i)+0.5_8*mesh%b, mesh%b)-0.5_8*mesh%b
2518 p(i) = p(i)/mesh%d + self%offset(i)
2519 if (
present(iv))
then 2520 p(i) = p(i) + mesh%h(iv)
2522 if (
present(roundup))
then 2523 ii(i) = ceiling(p(i))
2525 ii(i) = p(i) + 0.0001
2528 END FUNCTION index_only
2539 SUBROUTINE index_space2 (self, nbor, pos, iv, ii, p)
2541 class(
patch_t),
pointer :: nbor
2547 class(
mesh_t),
pointer :: m,mn
2550 real(8),
parameter:: eps=0.001_8
2551 integer,
save :: itimer=0
2552 real,
parameter:: sml = -1.0e-5
2559 pp = pp + pos(i) + m%h(iv)*(mn%d-m%d)
2560 if (self%periodic(i)) &
2561 pp = modulo(pp+0.5_8*m%b,m%b)-0.5_8*m%b
2565 ii(i)= floor(pp+eps)
2568 p(i) = pp -
real(ii(i),kind=4)
2569 if ((p(i)<0.0).and.(p(i)>sml)) p(i) = 0.0
2573 END SUBROUTINE index_space2
2578 FUNCTION index_only2 (self, nbor, pos, iv, roundup, closest)
RESULT (ii)
2580 class(
patch_t),
pointer :: nbor
2584 logical,
optional :: roundup
2585 logical,
optional :: closest
2587 class(
mesh_t),
pointer :: m, mn
2589 real(8),
parameter:: eps=0.001_8
2596 p(i) = p(i) + pos(i) + m%h(iv)*(mn%d-m%d)
2597 if (self%periodic(i)) &
2598 p(i) = modulo(p(i)+0.5_8*m%b,m%b)-0.5_8*m%b
2601 if (
present(roundup))
then 2602 ii(i) = ceiling(p(i)-eps)
2603 else if (
present(closest))
then 2606 ii(i) = floor(p(i)+eps)
2608 ii(i) = max(m%lb,min(m%ub,ii(i)))
2610 END FUNCTION index_only2
2615 FUNCTION count_edges (self, source)
RESULT (out)
2617 class(
patch_t),
pointer :: source
2621 dist = self%distance (source)
2624 if (abs(dist(i)) > self%ds(i)) out = out+1
2626 END FUNCTION count_edges
2633 SUBROUTINE interpolate (self, target, ii, iv, jt, pt)
2635 class(
patch_t),
pointer :: target
2636 integer :: ii(3), iv, jt(2)
2642 integer,
save :: itimer=0
2644 if (self%pervolume(iv))
then 2645 call self%interpolate_specific (
target, ii, iv, jt, pt)
2647 pos =
target%myposition (ii, iv)
2648 call self%index_space (pos, iv, j, p)
2654 if (self%unsigned(iv))
then 2655 target%mem(ii(1),ii(2),ii(3),iv,
target%it,1) = exp( &
2656 pt(1)*(q(3)*(q(2)*(q(1)*log(self%mem(j(1) ,j(2) ,j(3) ,iv,jt(1),1)) + &
2657 p(1)*log(self%mem(j(1)+1,j(2) ,j(3) ,iv,jt(1),1))) + &
2658 p(2)*(q(1)*log(self%mem(j(1) ,j(2)+1,j(3) ,iv,jt(1),1)) + &
2659 p(1)*log(self%mem(j(1)+1,j(2)+1,j(3) ,iv,jt(1),1)))) + &
2660 p(3)*(q(2)*(q(1)*log(self%mem(j(1) ,j(2) ,j(3)+1,iv,jt(1),1)) + &
2661 p(1)*log(self%mem(j(1)+1,j(2) ,j(3)+1,iv,jt(1),1))) + &
2662 p(2)*(q(1)*log(self%mem(j(1) ,j(2)+1,j(3)+1,iv,jt(1),1)) + &
2663 p(1)*log(self%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(1),1))))) + &
2664 pt(2)*(q(3)*(q(2)*(q(1)*log(self%mem(j(1) ,j(2) ,j(3) ,iv,jt(2),1)) + &
2665 p(1)*log(self%mem(j(1)+1,j(2) ,j(3) ,iv,jt(2),1))) + &
2666 p(2)*(q(1)*log(self%mem(j(1) ,j(2)+1,j(3) ,iv,jt(2),1)) + &
2667 p(1)*log(self%mem(j(1)+1,j(2)+1,j(3) ,iv,jt(2),1)))) + &
2668 p(3)*(q(2)*(q(1)*log(self%mem(j(1) ,j(2) ,j(3)+1,iv,jt(2),1)) + &
2669 p(1)*log(self%mem(j(1)+1,j(2) ,j(3)+1,iv,jt(2),1))) + &
2670 p(2)*(q(1)*log(self%mem(j(1) ,j(2)+1,j(3)+1,iv,jt(2),1)) + &
2671 p(1)*log(self%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(2),1))))))
2673 target%mem(ii(1),ii(2),ii(3),iv,
target%it,1) = &
2674 pt(1)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*self%mem(j(1) ,j(2) ,j(3) ,iv,jt(1),1) + &
2675 ( p(1))*self%mem(j(1)+1,j(2) ,j(3) ,iv,jt(1),1)) + &
2676 ( p(2))*((1.0-p(1))*self%mem(j(1) ,j(2)+1,j(3) ,iv,jt(1),1) + &
2677 ( p(1))*self%mem(j(1)+1,j(2)+1,j(3) ,iv,jt(1),1))) + &
2678 ( p(3))*((1.-p(2))*((1.0-p(1))*self%mem(j(1) ,j(2) ,j(3)+1,iv,jt(1),1) + &
2679 ( p(1))*self%mem(j(1)+1,j(2) ,j(3)+1,iv,jt(1),1)) + &
2680 ( p(2))*((1.0-p(1))*self%mem(j(1) ,j(2)+1,j(3)+1,iv,jt(1),1) + &
2681 ( p(1))*self%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(1),1)))) + &
2682 pt(2)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*self%mem(j(1) ,j(2) ,j(3) ,iv,jt(2),1) + &
2683 ( p(1))*self%mem(j(1)+1,j(2) ,j(3) ,iv,jt(2),1)) + &
2684 ( p(2))*((1.0-p(1))*self%mem(j(1) ,j(2)+1,j(3) ,iv,jt(2),1) + &
2685 ( p(1))*self%mem(j(1)+1,j(2)+1,j(3) ,iv,jt(2),1))) + &
2686 ( p(3))*((1.-p(2))*((1.0-p(1))*self%mem(j(1) ,j(2) ,j(3)+1,iv,jt(2),1) + &
2687 ( p(1))*self%mem(j(1)+1,j(2) ,j(3)+1,iv,jt(2),1)) + &
2688 ( p(2))*((1.0-p(1))*self%mem(j(1) ,j(2)+1,j(3)+1,iv,jt(2),1) + &
2689 ( p(1))*self%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(2),1))))
2692 if (
target%id==io%id_debug .and. io%verbose>2) &
2693 write(io%output,1)
target%id, ii, iv, self%id, j, jt, p, pt, &
2694 target%mem(ii(1),ii(2),ii(3),iv,
target%it,1)
2695 1
format(i6,1x,3i3,i4,2x,
"DBG patch_t%interpolate src:",i6, &
2696 " j,p;",3i4,1x,2i4,3f6.2,1x,2f6.2,
" out:",1p,e15.6)
2698 END SUBROUTINE interpolate
2705 SUBROUTINE interpolate_specific (self, target, ii, iv, jt, pt)
2707 class(
patch_t),
pointer :: target
2708 integer :: ii(3), iv, jt(2)
2710 real :: tmp(2,2,2,2)
2715 integer,
save :: itimer=0
2717 pos =
target%myposition (ii, iv)
2718 call self%index_space (pos, iv, j, p)
2723 associate(id=>self%idx%d)
2724 tmp(:,:,:,:) = self%mem(j(1):j(1)+1,j(2):j(2)+1,j(3):j(3)+1,iv,jt(1:2),1) &
2725 / self%mem(j(1):j(1)+1,j(2):j(2)+1,j(3):j(3)+1,id,jt(1:2),1)
2726 target%mem(ii(1),ii(2),ii(3),iv,
target%it,1) = &
2727 pt(1)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*tmp(1,1,1,1) + &
2728 ( p(1))*tmp(2,1,1,1)) + &
2729 ( p(2))*((1.0-p(1))*tmp(1,2,1,1) + &
2730 ( p(1))*tmp(2,2,1,1))) + &
2731 ( p(3))*((1.-p(2))*((1.0-p(1))*tmp(1,1,2,1) + &
2732 ( p(1))*tmp(2,1,2,1)) + &
2733 ( p(2))*((1.0-p(1))*tmp(1,2,2,1) + &
2734 ( p(1))*tmp(2,2,2,1)))) + &
2735 pt(2)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*tmp(1,1,1,2) + &
2736 ( p(1))*tmp(2,1,1,2)) + &
2737 ( p(2))*((1.0-p(1))*tmp(1,2,1,2) + &
2738 ( p(1))*tmp(2,2,1,2))) + &
2739 ( p(3))*((1.-p(2))*((1.0-p(1))*tmp(1,1,2,2) + &
2740 ( p(1))*tmp(2,1,2,2)) + &
2741 ( p(2))*((1.0-p(1))*tmp(1,2,2,2) + &
2742 ( p(1))*tmp(2,2,2,2))))
2746 target%mem(ii(1),ii(2),ii(3),iv,
target%it,1) = &
2747 target%mem(ii(1),ii(2),ii(3),iv,
target%it,1) * &
2748 target%mem(ii(1),ii(2),ii(3),id,
target%it,1)
2750 if (
target%id==io%id_debug .and. io%verbose>2) &
2751 write(io%output,1)
target%id, ii, iv, self%id, j, jt, p, pt, &
2752 target%mem(ii(1),ii(2),ii(3),iv,
target%it,1)
2753 1
format(i6,1x,3i3,i4,2x,
"DBG patch_t%interpolate_specific src:", &
2754 i6,
" j,p;",3i4,1x,2i4,3f6.2,1x,2f6.2,
" out:",1p,e15.6)
2756 END SUBROUTINE interpolate_specific
2758 SUBROUTINE update_position (self)
2761 if (all(self%velocity == 0.0))
return 2762 call trace%begin(
"patch_t%update_position")
2764 self%position = self%position + self%velocity * self%dtime
2765 self%llc_cart = self%llc_cart + self%velocity * self%dtime
2766 self%centre_cart = self%centre_cart + self%velocity * self%dtime
2768 self%mesh(:)%p = self%position(:)
2769 self%mesh(:)%llc_cart = self%llc_cart(:)
2772 where(self%position > 0.5*self%box)
2773 self%position = self%position - self%box
2774 self%mesh(:)%p = self%mesh(:)%p - self%box(:)
2775 self%llc_cart = self%llc_cart - self%box
2776 self%mesh(:)%llc_cart = self%mesh(:)%llc_cart - self%box
2777 self%centre_cart = self%centre_cart - self%box
2780 END SUBROUTINE update_position
2783 SUBROUTINE custom_refine (self)
2786 call io%abort (
"patch_t:: custom refine is triggered, but none provided. aborting.")
2787 END SUBROUTINE custom_refine
2792 FUNCTION get_overlap (self, patch, guards, l, u)
RESULT(overlap)
2794 class(
patch_t),
pointer:: patch
2795 integer:: l(3), u(3)
2796 logical:: guards, overlap
2798 real(8):: dist(3), la(3), ua(3)
2800 call trace%begin (
'download_t%overlap_range')
2801 overlap = self%overlaps (patch)
2807 dist = patch%position - self%position
2808 where (patch%periodic)
2809 dist = modulo(dist+0.5_8*patch%mesh%b,patch%mesh%b) - 0.5_8*patch%mesh%b
2817 la = dist - patch%size/2d0 - (patch%mesh%ng+1)*patch%mesh%d
2818 ua = dist + patch%size/2d0 + (patch%mesh%ng+1)*patch%mesh%d
2820 la = dist - patch%size/2d0
2821 ua = dist + patch%size/2d0
2823 l = floor(la/self%mesh%d + self%mesh%o)
2824 u = ceiling(ua/self%mesh%d + self%mesh%o)
2825 l = max(l,self%mesh%lb)
2826 u = min(u,self%mesh%ub)
2830 if (any(l > self%mesh%ui .or. u < self%mesh%li .or. l > u))
then 2831 write (stderr,
'(a,2(3x,a,i6,":",i2),3x,a,3f8.3,3x,a,2(3i4,2x))') &
2832 'patch_t%overlap_range false overlap, with ', &
2833 'self:', self%id, self%level, &
2834 'patch:', patch%id, patch%level, &
2835 'dist:', dist/self%ds, &
2842 where (self%mesh%n == 1)
2843 l = min(l,self%mesh%ub)
2844 u = max(u,self%mesh%lb)
2848 END FUNCTION get_overlap
2851 SUBROUTINE init_level (self)
2854 self%level = maxval(nint(log(self%box/self%ds) / &
2855 log(
real(self%refine_ratio))), mask=self%n > 1)
2856 END SUBROUTINE init_level
2859 Each thread uses a private timer data type, with arrays for start time and total time for each regist...
Support tic/toc timing, as in MATLAB, and accurate wallclock() function. The timing is generally much...
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
This module is a placeholder for shared information, solving the problem of making a piece of informa...
Boundary conditions for centered variables, which have the physical boundary in-between mesh points l...
Module holding anonymous pointers back to extras features.
Template module for mesh.
This index file has slot indices for all solver, all initially equal to zero It is the responsibility...
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...
Template module for tasks.