75 type(rt_boundaries_t) :: rt_boundaries
81 real,
pointer :: w_bin(:)
82 type(
rt_t),
pointer :: omega(:) => null()
83 type(
rt_t),
pointer :: main => null()
90 real(8) :: rt_llc(3)=0d0
91 real(8) :: rt_urc(3)=0d0
92 real(8) :: rt_timestep=0.01_8
93 real(8) :: warmup_end_time=0d0
94 integer :: n_mu=3, n_phi=4
96 logical :: vertical_propagation=.false.
97 logical :: warmup_done=.false.
98 real,
pointer :: q(:,:,:,:)
99 real,
pointer :: tt(:,:,:)
100 real,
pointer :: rk(:,:,:,:)
101 real,
pointer :: src(:,:,:,:)
102 real,
pointer :: qr(:,:,:)
105 procedure,
nopass:: init_task
106 procedure:: add_to_task_list
107 procedure:: init_omega
109 procedure:: is_upstream_of
110 procedure:: pre_update
111 procedure:: post_update
115 procedure:: setup_boundaries
116 procedure:: task_info
117 procedure:: nbor_info
118 procedure:: courant_condition
121 integer,
save:: verbose=0
122 type(
rt_t),
public:: rt
128 SUBROUTINE pre_init (self, link)
130 class(
link_t),
pointer :: link
132 integer:: dir, i_mu, i_phi, i_omega, n_omega, iostat, n
133 integer,
save:: n_mu=3, n_phi=4, n_bin=4, nt=3, n_warmup=3, ng=1
134 real,
save :: courant=0.5, cdtd=1.0
135 real,
save:: rt_grace=0.01
136 real(8),
save:: rt_llc(3)=-huge(1d0)
137 real(8),
save:: rt_urc(3)=+huge(1d0)
138 real(8),
save:: rt_timestep=0.01_8
139 logical,
save:: vertical_propagation=.false.
140 logical,
save:: first_time=.true.
141 logical,
save:: on=.true.
143 namelist /sc_rt_params/ on, verbose, nt, ng, rt_llc, rt_urc, rt_timestep, &
144 n_mu, n_phi, n_bin, rt_grace, vertical_propagation, n_warmup, courant, cdtd
148 call trace%begin (
'rt_t%pre_init')
153 read (io%input, sc_rt_params)
154 if (io%master)
write (io%output, sc_rt_params)
157 self%patch => gpatch%cast2gpatch(link%task)
161 self%n_warmup = n_warmup
163 self%rt_grace = rt_grace
166 self%rt_timestep = rt_timestep
167 self%vertical_propagation = vertical_propagation
170 self%courant = courant
172 call rt_integral%init
174 END SUBROUTINE pre_init
181 SUBROUTINE init_task (self)
182 class(
rt_t),
pointer:: self
185 class(
rt_t),
pointer :: omega
186 class(
link_t),
pointer :: link
187 class(*),
pointer :: dummy
188 integer:: n_omega, i_omega, dir, i_mu, i_phi, iostat, n, l(3), u(3), id_mpi
189 real,
dimension(:),
allocatable:: mu, w_mu
191 call trace%begin (
'rt_t%init_task')
207 n_omega = merge(1+self%n_phi*(self%n_mu-1), 0, self%n_mu > 0)
209 self%n_omega = n_omega
214 id_mpi = mpi%id%update(0)
215 self%id = id_mpi + (p%id-1)*(1+n_omega) + 1
217 write (io_unit%log,*)
'rt_t%init: id, mhd%id, mpi%id =', self%id, p%id, id_mpi
218 call self%patch_t%setup (position=p%position, size=p%size, n=p%n, nt=self%nt,&
219 ng=self%ng, nv=self%n_bin, nw=1)
226 self%origin = p%origin
227 self%llc_cart = p%llc_cart
228 self%llc_nat = p%llc_nat
229 self%position = p%position
230 self%centre_nat = p%position
231 self%restart = p%restart
232 self%warmup_end_time = p%time
236 self%t(self%it) = self%time
238 self%out_next = (int(self%time/io%out_time+1.0e-4)+1)*io%out_time
239 self%dt_fixed = self%rt_timestep
240 self%dtime = self%rt_timestep
241 self%dt = self%rt_timestep
242 self%grace = self%rt_grace
243 self%time = p%time - (self%n_warmup-self%grace)*self%dtime
244 self%eos_time = self%time - self%rt_timestep
245 self%iout = self%restart+1
246 self%periodic = p%periodic
247 self%verbose = verbose
248 call self%set (bits%no_density)
252 call self%setup_boundaries()
254 call self%rt_boundaries%init (self%mesh)
258 if (.not.
allocated(self%patch%heating_per_unit_volume))
then 259 allocate(self%patch%heating_per_unit_volume(p%gn(1), p%gn(2), p%gn(3)))
260 self%patch%heating_per_unit_volume = 0.0
265 allocate (self%w_bin (self%n_bin))
271 if (self%n_mu > 0)
then 272 allocate (mu(self%n_mu),w_mu(self%n_mu))
273 call radau%init (self%n_mu, mu, w_mu)
278 allocate (self%rk (self%gn(1),self%gn(2),self%gn(3),self%n_bin))
279 allocate (self%tt (self%gn(1),self%gn(2),self%gn(3)))
280 allocate (self%qr (self%gn(1),self%gn(2),self%gn(3)))
281 allocate (self%src(self%gn(1),self%gn(2),self%gn(3),self%n_bin))
282 call io%bits_mem(storage_size(self%tt ),product(shape(self%tt )),
'tt')
283 call io%bits_mem(storage_size(self%rk ),product(shape(self%rk )),
'rk')
284 call io%bits_mem(storage_size(self%qr ),product(shape(self%qr )),
'qr')
285 call io%bits_mem(storage_size(self%src ),product(shape(self%src)),
'src')
289 if (n_omega > 0)
allocate (self%omega(n_omega), source=self)
291 self%main => self%omega(1)%main
298 if (i_mu == self%n_mu) n = 1
302 omega => self%omega(i_omega)
303 omega%mu = dir*mu(i_mu)
304 omega%w = w_mu(i_mu)*math%pi4/n
305 omega%rt_phi = math%pi*(modulo(
real(2*(i_phi-1))/self%n_phi + &
306 (dir+1)/2.0, 2.0) - 1.0)
307 if (i_mu == self%n_mu) omega%rt_phi = 0.0
308 call self%init_omega(i_omega)
315 if (verbose > 0 .and. self%id == 1)
then 317 omega => self%omega(i_omega)
319 write(io%output,
'(a,3i7,f8.3)') &
320 'rt_t%init: i_omega, elevation, phi, mu =', i_omega, &
321 nint(asin(omega%mu)*180./math%pi), nint(omega%rt_phi*180./math%pi), &
329 l = p%mesh%li - self%mesh%ng
330 u = p%mesh%ui + self%mesh%ng
331 self%mem => p%mem(l(1):u(1),l(2):u(2),l(3):u(3),:,:,:)
340 END SUBROUTINE init_task
345 SUBROUTINE add_to_task_list (self)
347 class(
list_t),
pointer:: task_list
351 call trace%begin (
'rt_t%add_tasks_to')
353 task_list => self%patch%task_list
354 if (self%is_clear(bits%no_rt))
then 355 call task_list%prepend_link (self%link)
356 do i_omega=1,self%n_omega
357 call task_list%prepend_link (self%omega(i_omega)%link)
361 END SUBROUTINE add_to_task_list
366 SUBROUTINE init_omega (self, i_omega)
367 class(
rt_t),
target:: self
370 class(
rt_t),
pointer:: omega
374 call trace%begin (
'rt_t%init_omega')
375 omega => self%omega(i_omega)
378 omega%id = self%id + i_omega
379 call omega%task_t%init
381 write (io_unit%log,*)
'rt_t%init_omega: id, self%id, mhd%id =', omega%id, &
382 self%id, self%patch%id
384 omega%i_omega = i_omega
385 omega%it = max(omega%nt-1,1)
388 omega%nv = self%n_bin
390 allocate (omega%unsigned(omega%nv), omega%pervolume(omega%nv))
391 omega%unsigned = .false.
392 omega%pervolume = .false.
394 allocate (omega%t(omega%nt), omega%dt(omega%nt), omega%iit(omega%nt))
399 omega%time = self%time + 2.0*self%grace*self%dtime
401 omega%dt = omega%dtime
403 allocate (omega%mem(self%gn(1),self%gn(2),self%gn(3),self%n_bin,self%nt,1))
404 call io%bits_mem(storage_size(omega%mem) ,product(shape(omega%mem)) ,
'q')
405 omega%q => omega%mem(:,:,:,:,omega%it,1)
410 omega%rt_hat(1) = sin(acos(omega%mu))*cos(omega%rt_phi)
411 omega%rt_hat(2) = sin(acos(omega%mu))*sin(omega%rt_phi)
412 omega%rt_hat(3) = omega%mu
413 loc = maxloc(abs(omega%rt_hat))
418 allocate (omega%link)
419 omega%link%task => omega
421 END SUBROUTINE init_omega
428 FUNCTION needs (self, omega, task)
430 class(
rt_t):: self, omega
433 real:: ds1(3), ds2(3), sinth
434 real,
parameter:: eps=1e-4
435 integer:: is1(3), is2(3), i, lbit(3), ubit(3)
436 logical:: lower, upper, upstream
441 ds1 = self%position-task%position
442 where (self%periodic)
443 ds1 = modulo(ds1 + 0.5_8*self%mesh%b, self%mesh%b) - 0.5_8*self%mesh%b
446 sinth = merge(0.0, sqrt(1.0-omega%mu**2), abs(omega%mu)==1.0)
447 ds2 = [sinth*cos(omega%rt_phi),sinth*sin(omega%rt_phi),omega%mu]
451 if (ds1(i) > +eps) is1(i)=1
452 if (ds2(i) > +eps) is2(i)=1
453 if (ds1(i) < -eps) is1(i)=-1
454 if (ds2(i) < -eps) is2(i)=-1
460 if (task%is_set(bits%no_rt))
then 467 lbit = [bits%xl, bits%yl, bits%zl]
468 ubit = [bits%xu, bits%yu, bits%zu]
470 lower = self%boundaries%is_set(lbit(i))
471 upper = self%boundaries%is_set(ubit(i))
472 if (upper .and. ds1(i)==1)
then 475 if (lower .and. ds1(i)==-1)
then 482 if (abs(omega%mu)==1.0)
then 483 needs = all(is1==is2)
487 else if (any(is1*is2 < 0))
then 492 else if (is2(1)==0 .and. is1(1)/=0)
then 494 else if (is2(2)==0 .and. is1(2)/=0)
then 497 upstream = task%is_upstream_of(omega)
498 if (needs .neqv. upstream)
then 499 print
'(a,2i4,2f8.3,l4,2(3f8.3,2x),2i4)', &
500 'needs and is_upstream of disagree on', omega%id, task%id, omega%mu, &
501 omega%rt_phi, needs,
real(omega%task_t%distance(task),kind=4), &
502 omega%mesh%s, omega%i_omega, task%i_omega
519 LOGICAL FUNCTION is_upstream_of (self, other)
520 class(
rt_t):: self, other
521 real(8):: p(3), d(3), dh(3)
528 d = other%distance (self)
530 is_upstream_of = all(dh > -self%ds)
535 print *, self%id, other%id
536 print 1,
'sp:', self%position ,
'sm:', self%mesh%p
537 print 1,
'op:', other%position ,
'om:', other%mesh%p
539 print 1,
'dh:', dh ,
'rt:', self%rt_hat, is_upstream_of
540 1
format(2(a,3f8.3,2x),l3)
545 if (is_upstream_of)
then 548 associate(r => self%mesh(i)%r)
553 if (abs(dh(i)) > self%ds(i))
then 554 if (d(i) > +self%ds(i))
then 555 p(i) = p(i) + r(self%mesh(i)%uo)
556 else if (d(i) < -self%ds(i))
then 557 p(i) = p(i) + r(self%mesh(i)%lo)
565 if (.not. other%overlaps_point (p))
then 566 is_upstream_of = .false.
570 print
'(3i4,l4,5(2x,a,3f8.3))', self%id, other%id, self%axis, is_upstream_of, &
571 'sp:',self%mesh%p,
'op:',other%mesh%p,
'd:',abs(d), &
572 's2:',other%mesh%s/2d0
573 END FUNCTION is_upstream_of
580 SUBROUTINE setup_boundaries(self)
583 real(8) :: llc(3), urc(3)
584 integer :: lbit(3), ubit(3), i
586 call self%init_bdries
587 llc = self%position-0.5_8*self%size
588 urc = self%position+0.5_8*self%size
590 lbit = [bits%xl, bits%yl, bits%zl]
591 ubit = [bits%xu, bits%yu, bits%zu]
594 if ((urc(i) < self%rt_llc(i)))
then 595 call self%patch%set(bits%no_rt)
597 else if (llc(i) <= self%rt_llc(i))
then 598 call self%boundaries%set (lbit(i))
599 self%mesh(i)%lower_boundary = .true.
602 if ((llc(i) > self%rt_urc(i)))
then 603 call self%patch%set(bits%no_rt)
605 else if (urc(i) >= self%rt_urc(i))
then 606 call self%boundaries%set (ubit(i))
607 self%mesh(i)%upper_boundary = .true.
610 END SUBROUTINE setup_boundaries
617 SUBROUTINE pre_update (self)
620 class(
rt_t),
pointer:: omega
621 integer:: l(3), u(3), i_omega, i_bin
622 integer,
save:: itimer=0
626 call trace%begin (
'rt_t%pre_update', itimer=itimer)
631 do i_omega=1,self%n_omega
632 omega => self%omega(i_omega)
633 do i_bin=1,self%n_bin
634 self%qr = self%qr + self%rk(:,:,:,i_bin)*self%w_bin(i_bin)*omega%w* &
635 omega%mem(:,:,:,i_bin,omega%it,1)
638 call self%patch%aux%register (
'qr', self%qr)
640 l = self%patch%mesh%li - self%mesh%ng
641 u = self%patch%mesh%ui + self%mesh%ng
642 self%patch%heating_per_unit_volume(l(1):u(1),l(2):u(2),l(3):u(3)) = &
643 self%patch%heating_per_unit_volume(l(1):u(1),l(2):u(2),l(3):u(3)) + self%qr
645 call trace%end (itimer)
646 END SUBROUTINE pre_update
653 SUBROUTINE post_update (self)
656 class(gpatch_t),
pointer:: mhd
657 class(
rt_t),
pointer:: omega
661 if (verbose == -1)
then 663 timer%n_mhd = timer%n_mhd+1
664 write (io_unit%queue,*)
'mhd', timer%n_mhd, wallclock(), &
665 self%patch%id, self%id, mhd%time+mhd%dtime, self%on, self%rotated
666 if (timer%n_mhd > 1000)
then 673 if (.not.self%on)
return 674 call trace%begin (
'rt_t%post_update')
679 self%dtime = mhd%dtime
680 self%time = mhd%time + self%grace*self%dtime
682 omega => self%omega(i)
683 omega%time = mhd%time
689 if (self%vertical_propagation)
then 690 omega%dtime = mhd%dtime*abs(omega%mu)
692 omega%dtime = mhd%dtime + 2*self%grace*self%dtime
697 END SUBROUTINE post_update
705 SUBROUTINE update (self)
708 class(
rt_t),
pointer :: omega
709 integer:: i_omega, i_bin
710 real,
dimension(:,:,:),
pointer:: q
712 integer,
save:: itimer=0
716 call trace%begin (
'rt_t%update', itimer=itimer)
721 if (self%i_omega==0)
then 722 call self%rt_eos (self%time+self%dtime)
728 self%q => self%mem(:,:,:,:,self%new,1)
730 if (verbose == -1)
then 732 timer%n_solve = timer%n_solve+1
733 write (io_unit%queue,*)
'ome', timer%n_solve, wallclock(), &
734 self%patch%id, self%id, self%time+self%dtime, self%i_omega
737 call trace%end (itimer)
738 END SUBROUTINE update
744 SUBROUTINE rt_eos (self, time)
749 class(gpatch_t),
pointer :: patch
750 integer,
save:: itimer=0
751 real,
dimension(:,:,:),
pointer:: d, e, pgas, fmax
753 if (time > self%eos_time)
then 755 if (time > self%eos_time)
then 759 write (io_unit%queue,*)
'eos', timer%n_mhd, wallclock(), &
760 self%patch%id, self%id, time
761 else if (verbose == -1)
then 762 write (io_unit%queue,*)
'eo0', timer%n_mhd, wallclock(), &
763 self%patch%id, self%id, time, self%eos_time
775 call trace%begin (
'rt_t%rt_eos', itimer=itimer)
778 allocate (pgas(m(1),m(2),m(3)),fmax(m(1),m(2),m(3)))
780 d => self%mem(:,:,:,patch%idx%d,patch%it,1)
781 e => self%mem(:,:,:,patch%idx%e,patch%it,1)
782 if (verbose > 2)
then 783 call self%stats (d,
'd')
784 call self%stats (e,
'e')
786 call eos%lookup_table (m, d=d, e=e, rk=self%rk, src=self%src, tt=self%tt, pg=pgas)
791 stefan = cgs%stefan/(scaling%p*scaling%u)*scaling%temp**4
792 do ib=1, eos%n_lambda
793 fmax = math%pi2*self%src(:,:,:,ib)/pgas*(self%rk(:,:,:,ib)*minval(self%ds))
795 fmax = fmax/(1.+(self%rk(:,:,:,ib)*minval(self%ds))**2)
796 u_max = 16./3.*self%fmaxval(fmax)*self%courant/self%cdtd*self%epsilon
797 self%u_max = max(self%u_max,u_max)
799 self%omega(:)%u_max = self%u_max
800 if (verbose > 2)
then 801 call self%stats (self%tt ,
'tt')
802 call self%stats (pgas,
'pgas')
803 call self%stats (self%rk(:,:,:,1),
'rkap')
804 call self%stats (fmax,
'fmax')
807 deallocate (fmax, pgas)
808 if (.not.self%warmup_done)
then 809 if (self%time > self%warmup_end_time)
then 810 self%warmup_done=.true.
814 call trace%end (itimer)
815 end subroutine omp_eos
816 END SUBROUTINE rt_eos
821 SUBROUTINE solve (self)
824 integer,
save:: itimer=0
826 if (self%i_omega==0)
return 827 call trace%begin (
'rt_t%solve', itimer=itimer)
829 call self%rt_boundaries%condition (mem=self%q, src=self%src, rk=self%rk, &
830 tt=self%tt, mu=self%rt_hat)
831 call rt_integral%solve (self%rk, self%src, self%q, self%mesh, self%axis, &
833 call self%patch%aux%register (
'rkap', self%rk )
834 call self%patch%aux%register (
'src', self%src)
835 call trace%end (itimer)
841 SUBROUTINE task_info (self)
844 write(io_unit%mpi,
'(1x,a,12x,i6,2i4,7x,2l1,2x,3f9.3,i4,2f7.3)') &
845 'rt_t%task_info: id, rank, level, BV, pos, i_omega, mu, phi =', &
846 self%id, self%rank, self%level, &
847 self%is_set(bits%boundary), self%is_set(bits%virtual), self%position, &
848 self%i_omega, self%mu, self%rt_phi
849 END SUBROUTINE task_info
854 SUBROUTINE nbor_info (self, nbor)
858 write(io_unit%mpi,
'(3x,a,i6,2i4,2x,3l1,2x,2l1,2x,3f9.3,i4,2f7.3)') &
859 'rt_t%nbor_info: id, rank, level, needed, needs_me, download, BV, pos =', &
860 self%id, self%rank, self%level, nbor%needed, nbor%needs_me, nbor%download, &
861 self%is_set(bits%boundary), self%is_set(bits%virtual), &
862 self%position, self%i_omega, self%mu, self%rt_phi
863 END SUBROUTINE nbor_info
868 SUBROUTINE courant_condition (self, detailed_timer)
870 logical,
optional:: detailed_timer
872 call trace%begin (
'rt_t%courant_condition')
873 self%patch%u_max = max(self%patch%u_max, self%u_max)
875 END SUBROUTINE courant_condition
880 SUBROUTINE output (self)
883 END SUBROUTINE output
Each thread uses a private timer data type, with arrays for start time and total time for each regist...
Given a direction Omega, which may align mainly with the x-. y-, or z-axis, solve the integral form o...
The gpath_t layer now essentially only handles restarts.
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.
Module to set up angles and bins for an RT patch with only absorption. This means that the only extra...
Equation of state module for any sort of tables, provided by a reader.
Fundamental constants in CGS and SI units.
Boundary conditions for centered variables, which have the physical boundary in-between mesh points l...
Fundamental constants in CGS and SI units.
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.
Module with which one can register any number of pointers to real or integer arrays, and then output the contents of the links later.
Template module for tasks.