DISPATCH
rt_mod.f90
1 !===============================================================================
2 !> Module to set up angles and bins for an RT patch with only absorption.
3 !> This means that the only extra memory with time slot that is needed is
4 !> one slot for the net radiation heating, \Sum_bins \Sum_angles Q_bin_angle
5 !>
6 !> The data structure consists of self, an extension of the gpatch_t data
7 !> type, which also contains n_omega clones of itself, with the specific info
8 !> for each ray direction. Since the data type is an extension of patch_t, we
9 !> can use the standard guard zone loading procedures in download_mod or
10 !> data_hub_mod to fill the RT guard zones, and since the sub-tasks and main
11 !> task have the same data type, the sub-tasks can be cloned from the main task.
12 !>
13 !> Since directional information becomes available at different times, the
14 !> guard zone filling must take place for each angle, just before the solver
15 !> is called.
16 !>
17 !> The self%time, self%rt_timestep, and self%rt_next control the scheduling of
18 !> the tasks. The MHD tasks require that the times of their nbors are "ahead",
19 !> allowing for a grace fraction of the nbor time step. The RT sub-tasks
20 !> advance their time by self%rt_timestep (effectively their self%dtime) when
21 !> done, hence they must update when their MHD master task passes their own time,
22 !> which conforms with standard nbor relations. To clarify, with the times
23 !>
24 !> MHD time: mhd%time
25 !> RT time: mhd%rt%time
26 !> ray time: mhd%rt%omega%time
27 !>
28 !> The RT_time is a "rabbit time", which gives the time until which the MHD task
29 !> should still be happy with the RT solution. When the parent MHD task passes
30 !> that time, the RT sub-tasks become out-of-date, and they perform solutions.
31 !> Afterwards, they advance their "rabbit time" again, and when all of them have
32 !> done so, the main RT task becomes ready to update, and when it runs it sums
33 !> up the sub-task contributions, and advances its own "rabbit time", which then
34 !> allows the MHD task to update. Clearly, the EOS call that gives values to
35 !> the RT tasks at the new MHD time has to be made after the MHD task has passed
36 !> by, but before the first sub-task executes. This could be detected either in
37 !> the MHD task, or in the RT sub-tasks, where the first detection should lead
38 !> to locking the EOS while executing it, followed by another test (so a double
39 !> critical region type of thing).
40 !>
41 !> Initially, we wish that only the RT states are ready to update, and we set
42 !> things up to compute the EOS, solve the RT, sum up Q, and only then let the
43 !> MHD update. This means that initially, the EOS_time should be negative, the
44 !> RT times should be RT_timestep/2 negative, and the MHD_time should be zero.
45 !===============================================================================
46 MODULE rt_mod
47  USE rt_boundaries_mod
48  USE rt_integral_mod
49  USE radau_mod
50  USE boundaries_mod
51  USE eos_mod
52  USE list_mod
53  USE gpatch_mod
54  USE link_mod
55  USE task_mod
56  USE index_mod
57  USE mesh_mod
58  USE io_mod
59  USE aux_mod
60  USE timer_mod
61  USE scaling_mod
62  USE trace_mod
63  USE mpi_mod
64  USE omp_mod
65  USE omp_timer_mod
66  USE math_mod
67  USE bits_mod
68  implicit none
69  private
70  !-----------------------------------------------------------------------------
71  ! RT patch data type, extended from gpatch, with extra RT memory and parameters
72  !-----------------------------------------------------------------------------
73  type, public, extends(gpatch_t):: rt_t
74  class(gpatch_t), pointer :: patch ! pointer to MHD patch
75  type(rt_boundaries_t) :: rt_boundaries ! RT boundary conditions
76  integer :: n_warmup=5 ! number of initial solutions
77  integer :: n_bin=4 ! number of frequency bins
78  integer :: n_omega=0 ! number of space-angles
79  integer :: i_omega=0 ! space-angle index
80  integer :: axis=0 ! neares axis
81  real, pointer :: w_bin(:) ! frequency weight
82  type(rt_t), pointer :: omega(:) => null() ! different RT directions
83  type(rt_t), pointer :: main => null() ! point to main RT task
84  real :: cdtd = 1.0
85  real :: mu, rt_phi, w ! mu, phi, weights
86  real :: rt_hat(3) ! RT direction
87  real :: epsilon =1.0 ! scattering coefficient
88  real(8) :: eos_time ! laste EOS time
89  real :: rt_grace=0.01 ! deal with roundoff
90  real(8) :: rt_llc(3)=0d0 ! lower left corner
91  real(8) :: rt_urc(3)=0d0 ! upper right corner
92  real(8) :: rt_timestep=0.01_8 ! constant time step
93  real(8) :: warmup_end_time=0d0 ! time when RT time should start varying
94  integer :: n_mu=3, n_phi=4 ! no. of angles
95  logical :: on=.false. ! allows local turn-off
96  logical :: vertical_propagation=.false. ! see code
97  logical :: warmup_done=.false.
98  real, pointer :: q(:,:,:,:) ! net radiative heating
99  real, pointer :: tt(:,:,:) ! EOS temperature
100  real, pointer :: rk(:,:,:,:) ! EOS rho*kappa
101  real, pointer :: src(:,:,:,:) ! EOS source function
102  real, pointer :: qr(:,:,:) ! net radiative heating
103  contains
104  procedure:: pre_init ! read parameters
105  procedure, nopass:: init_task ! init RT task
106  procedure:: add_to_task_list ! add to task list
107  procedure:: init_omega ! init omega task
108  procedure:: needs ! check upstream
109  procedure:: is_upstream_of ! check upstream
110  procedure:: pre_update ! pre MHD
111  procedure:: post_update ! post MHD
112  procedure:: update ! RT task update
113  procedure:: rt_eos ! EOS for RT
114  procedure:: solve ! RT solver
115  procedure:: setup_boundaries ! BCs
116  procedure:: task_info ! formatted info print
117  procedure:: nbor_info ! formatted info print
118  procedure:: courant_condition ! amend Courant condition
119  procedure:: output ! intercept output
120  end type
121  integer, save:: verbose=0 ! verbosity
122  type(rt_t), public:: rt ! public instance
123 CONTAINS
124 
125 !===============================================================================
126 !> Read parameters
127 !===============================================================================
128 SUBROUTINE pre_init (self, link)
129  class(rt_t):: self
130  class(link_t), pointer :: link
131  !-----------------------------------------------------------------------------
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.
142  !-----------------------------------------------------------------------------
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
145  !-----------------------------------------------------------------------------
146  ! Read parameters
147  !-----------------------------------------------------------------------------
148  call trace%begin ('rt_t%pre_init')
149  !$omp critical (input_cr)
150  if (first_time) then
151  first_time = .false.
152  rewind(io%input)
153  read (io%input, sc_rt_params)
154  if (io%master) write (io%output, sc_rt_params)
155  end if
156  !$omp end critical (input_cr)
157  self%patch => gpatch%cast2gpatch(link%task)
158  self%on = on
159  self%nt = nt
160  self%ng = ng
161  self%n_warmup = n_warmup
162  self%n_bin = n_bin
163  self%rt_grace = rt_grace
164  self%rt_llc = rt_llc
165  self%rt_urc = rt_urc
166  self%rt_timestep = rt_timestep
167  self%vertical_propagation = vertical_propagation
168  self%n_mu = n_mu
169  self%n_phi = n_phi
170  self%courant = courant
171  self%cdtd = cdtd
172  call rt_integral%init
173  call trace%end()
174 END SUBROUTINE pre_init
175 
176 !===============================================================================
177 !> Initialize the rt_t data types. To setup for n_warmup initial solutions, we
178 !> set the time rt_t times = -n_warmup * rt_t%timestep, making sure to call
179 !> rt_t%rt_eos for the initial MHD state.
180 !===============================================================================
181 SUBROUTINE init_task (self)
182  class(rt_t), pointer:: self
183  ! ............................................................................
184  class(gpatch_t), pointer :: p
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
190  ! ----------------------------------------------------------------------------
191  call trace%begin ('rt_t%init_task')
192  p => self%patch
193  !-----------------------------------------------------------------------------
194  ! interconnect itself
195  !-----------------------------------------------------------------------------
196  nullify(self%link)
197  allocate(link)
198  link%task => self
199  self%link => link
200  self%track = .false.
201  self%type = 'rt_t'
202  self%kind = 'rt_sc'
203  self%box = p%box
204  !-----------------------------------------------------------------------------
205  ! Number of ray directions
206  !-----------------------------------------------------------------------------
207  n_omega = merge(1+self%n_phi*(self%n_mu-1), 0, self%n_mu > 0)
208  n_omega = 2*n_omega
209  self%n_omega = n_omega
210  !-----------------------------------------------------------------------------
211  ! Set the RT task ID in a predictable way, spacing by (1+n_omega) behind the
212  ! largest MPI-wide ID used so far, to make room for the omega tasks
213  !-----------------------------------------------------------------------------
214  id_mpi = mpi%id%update(0)
215  self%id = id_mpi + (p%id-1)*(1+n_omega) + 1
216  if (verbose == -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)
220  !-----------------------------------------------------------------------------
221  ! Synchronize with MHD patch
222  ! ----------------------------------------------------------------------------
223  self%rank = p%rank
224  self%level = p%level
225  self%size = p%size
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
233  self%iit = 1
234  self%mu = 1.0
235  self%rt_phi = 0.0
236  self%t(self%it) = self%time
237  ! -- add an `eps` to avoid round-off errors
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)
249  !-----------------------------------------------------------------------------
250  ! Check if this MHD patch needs RT
251  !-----------------------------------------------------------------------------
252  call self%setup_boundaries()
253  if (self%on) then
254  call self%rt_boundaries%init (self%mesh)
255  !---------------------------------------------------------------------------
256  ! Allocate heating for MHD
257  !---------------------------------------------------------------------------
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
261  end if
262  !---------------------------------------------------------------------------
263  ! Integration weights (w_bin is baked into the source function in this case)
264  !---------------------------------------------------------------------------
265  allocate (self%w_bin (self%n_bin))
266  self%w_bin = 1.0
267  !---------------------------------------------------------------------------
268  ! Assign angles and integration weight and allocate the required RT solvers
269  ! for all angles
270  !---------------------------------------------------------------------------
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)
274  end if
275  !---------------------------------------------------------------------------
276  ! Allocate data, belonging to the main task
277  !---------------------------------------------------------------------------
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')
286  !---------------------------------------------------------------------------
287  ! Allocate n_omega copies of self, cloning self
288  !---------------------------------------------------------------------------
289  if (n_omega > 0) allocate (self%omega(n_omega), source=self)
290  self%i_omega = 0
291  self%main => self%omega(1)%main
292  !---------------------------------------------------------------------------
293  ! Assign angles and initialize clones
294  !----------------------------------------------------------------------------
295  i_omega = 0
296  do i_mu=1,self%n_mu
297  n = self%n_phi
298  if (i_mu == self%n_mu) n = 1
299  do i_phi=1,n
300  do dir=-1,1,2
301  i_omega = i_omega+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)
309  end do
310  end do
311  end do
312  !---------------------------------------------------------------------------
313  ! Check task angles
314  !---------------------------------------------------------------------------
315  if (verbose > 0 .and. self%id == 1) then
316  do i_omega=1,n_omega
317  omega => self%omega(i_omega)
318  if (verbose > 0) &
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), &
322  omega%mu
323  end do
324  end if
325  !---------------------------------------------------------------------------
326  ! Point mem to a patch mem section -- download_mod should be able handle a
327  ! target with fewer guard cells than the source
328  !---------------------------------------------------------------------------
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),:,:,:)
332  else
333  !---------------------------------------------------------------------------
334  ! Deallocate the self%mem for this task, but keep the task, for simplicity
335  !---------------------------------------------------------------------------
336  call self%dealloc()
337  self%n_omega = 0
338  end if
339  call trace%end()
340 END SUBROUTINE init_task
341 
342 !===============================================================================
343 !> Add the (1+n_omega) RT tasks to the task list
344 !===============================================================================
345 SUBROUTINE add_to_task_list (self)
346  class(rt_t):: self
347  class(list_t), pointer:: task_list
348  !.............................................................................
349  integer:: i_omega
350  !-----------------------------------------------------------------------------
351  call trace%begin ('rt_t%add_tasks_to')
352  !-----------------------------------------------------------------------------
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)
358  end do
359  end if
360  call trace%end()
361 END SUBROUTINE add_to_task_list
362 
363 !===============================================================================
364 !> Initialize rt sub-tasks; note that all self values are initially cloned
365 !===============================================================================
366 SUBROUTINE init_omega (self, i_omega)
367  class(rt_t), target:: self
368  integer:: i_omega
369  ! ............................................................................
370  class(rt_t), pointer:: omega
371  real:: mu(3)
372  integer:: loc(1)
373  ! ----------------------------------------------------------------------------
374  call trace%begin ('rt_t%init_omega')
375  omega => self%omega(i_omega)
376  omega%main => self
377  ! -- assign unique and globally valid ID
378  omega%id = self%id + i_omega
379  call omega%task_t%init
380  if (verbose == -1) &
381  write (io_unit%log,*) 'rt_t%init_omega: id, self%id, mhd%id =', omega%id, &
382  self%id, self%patch%id
383  ! -- set values that should be different from self
384  omega%i_omega = i_omega
385  omega%it = max(omega%nt-1,1)
386  omega%new = omega%nt
387  omega%nw = 1
388  omega%nv = self%n_bin
389  ! -- allocate unsigned and pervolume
390  allocate (omega%unsigned(omega%nv), omega%pervolume(omega%nv))
391  omega%unsigned = .false.
392  omega%pervolume = .false.
393  ! -- allocate time slots
394  allocate (omega%t(omega%nt), omega%dt(omega%nt), omega%iit(omega%nt))
395  !-----------------------------------------------------------------------------
396  ! Setting omega%time beyond self%time to trigger first call to self%update(),
397  ! which is needed to compute the EOS data
398  !-----------------------------------------------------------------------------
399  omega%time = self%time + 2.0*self%grace*self%dtime
400  omega%t = omega%time
401  omega%dt = omega%dtime
402  omega%iit = 1
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)
406  omega%mem = 0.0
407  !-----------------------------------------------------------------------------
408  ! Determine closest axis direction
409  !-----------------------------------------------------------------------------
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))
414  omega%axis = loc(1)
415  !-----------------------------------------------------------------------------
416  ! Allocate task link and set link%task
417  !-----------------------------------------------------------------------------
418  allocate (omega%link)
419  omega%link%task => omega
420  call trace%end()
421 END SUBROUTINE init_omega
422 
423 !===============================================================================
424 !> Determine if a task is needed (is upstream). We are being presented with
425 !> the nbors of the MHD task, so there is no task above the top-most task, but
426 !> there may well be tasks below the rt_bottom, which are not needed as RT nbors
427 !===============================================================================
428 FUNCTION needs (self, omega, task)
429  logical:: needs
430  class(rt_t):: self, omega
431  class(rt_t):: task
432  !.............................................................................
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
437  !-----------------------------------------------------------------------------
438  ! If self is at a boundary and omega points inward, or the task has no RT,
439  ! then the task is outside and is not needed
440  !-----------------------------------------------------------------------------
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
444  end where
445  ds1 = ds1/self%box
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]
448  is1 = 0
449  is2 = 0
450  do i=1,3
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
455  end do
456  needs = .true.
457  !-----------------------------------------------------------------------------
458  ! If the task does not have RT it is not needed
459  !-----------------------------------------------------------------------------
460  if (task%is_set(bits%no_rt)) then
461  needs = .false.
462  return
463  end if
464  !-----------------------------------------------------------------------------
465  ! If the nbor lies outside the boundaries it is not needed
466  !-----------------------------------------------------------------------------
467  lbit = [bits%xl, bits%yl, bits%zl]
468  ubit = [bits%xu, bits%yu, bits%zu]
469  do i=1,3
470  lower = self%boundaries%is_set(lbit(i))
471  upper = self%boundaries%is_set(ubit(i))
472  if (upper .and. ds1(i)==1) then
473  needs = .false.
474  end if
475  if (lower .and. ds1(i)==-1) then
476  needs = .false.
477  end if
478  end do
479  !-----------------------------------------------------------------------------
480  ! Vertical rays depend only on exact alignment
481  !-----------------------------------------------------------------------------
482  if (abs(omega%mu)==1.0) then
483  needs = all(is1==is2)
484  !-----------------------------------------------------------------------------
485  ! For inclined rays, define upstream as "not downstream in any direction"
486  !-----------------------------------------------------------------------------
487  else if (any(is1*is2 < 0)) then
488  needs = .false.
489  !-----------------------------------------------------------------------------
490  ! If a horizontal ray component is axis aligned, exclude non-aligned nbors
491  !-----------------------------------------------------------------------------
492  else if (is2(1)==0 .and. is1(1)/=0) then
493  needs = .false.
494  else if (is2(2)==0 .and. is1(2)/=0) then
495  needs = .false.
496  end if
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
503  !needs = upstream
504 end if
505 END FUNCTION needs
506 
507 !===============================================================================
508 !> Determine if the self task is "upstream" of the other task; i.e., if its
509 !> guard cells that lie between the two patches overlap with the internal
510 !> region of the other patch
511 !>
512 !> In the most general case one needs to first identify the "faces" of the task
513 !> that are at play. Those correspond to the components of the ray direction
514 !> that are non-zero. For each upstream candidate, one should run through the
515 !> candidate faces, identified by squares (or more generally flat volumes), and
516 !> check if they have overlap in all three directions with the internal region
517 !> of the candidate upstream task.
518 !===============================================================================
519 LOGICAL FUNCTION is_upstream_of (self, other)
520  class(rt_t):: self, other
521  real(8):: p(3), d(3), dh(3)
522  integer:: i
523  logical:: debug
524  !-----------------------------------------------------------------------------
525  ! Check if the direction vector d is in the same octant; to be correct for
526  ! periodic meshes this must use the patch_t%distance() function!
527  !-----------------------------------------------------------------------------
528  d = other%distance (self) ! direction towards other
529  dh = d*self%rt_hat ! projection onto rt_hat
530  is_upstream_of = all(dh > -self%ds) ! all components "not negative"
531 debug=.false.
532 !debug = self%position(2)==other%position(2) !&
533  !.and. abs(self%position(2)-0.5d0) < 0.01
534 if (debug) then
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
538 print 1, ' d:', d
539 print 1, 'dh:', dh , 'rt:', self%rt_hat, is_upstream_of
540 1 format(2(a,3f8.3,2x),l3)
541 end if
542  !-----------------------------------------------------------------------------
543  ! If we have a possible upstream task
544  !-----------------------------------------------------------------------------
545  if (is_upstream_of) then
546  p = self%position
547  do i=1,3
548  associate(r => self%mesh(i)%r)
549  !-------------------------------------------------------------------------
550  ! If the direction vector projected on the ray is non-zero (w roundoff),
551  ! then move the point to the guard zone in the other task direction
552  !-------------------------------------------------------------------------
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)
558  end if
559  end if
560  end associate
561  end do
562  !---------------------------------------------------------------------------
563  ! If that doesn't overlap with the other task, self is not upstream
564  !---------------------------------------------------------------------------
565  if (.not. other%overlaps_point (p)) then
566  is_upstream_of = .false.
567  end if
568  end if
569 if (debug) &
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
574 
575 !===============================================================================
576 !> Check if patch belongs to an RT box or is at its boundary. The mesh flags
577 !> lower_boundary and upper_boundary are set for the tasks that contain bndries,
578 !> and self%on is set false for tasks that are entirely outside.
579 !===============================================================================
580 SUBROUTINE setup_boundaries(self)
581  class(rt_t):: self
582  ! ----------------------------------------------------------------------------
583  real(8) :: llc(3), urc(3)
584  integer :: lbit(3), ubit(3), i
585  ! ............................................................................
586  call self%init_bdries
587  llc = self%position-0.5_8*self%size
588  urc = self%position+0.5_8*self%size
589  !
590  lbit = [bits%xl, bits%yl, bits%zl] ! lower boundary bits -- convenience
591  ubit = [bits%xu, bits%yu, bits%zu] ! upper boundary bits -- convenience
592  !
593  do i = 1,3
594  if ((urc(i) < self%rt_llc(i))) then
595  call self%patch%set(bits%no_rt)
596  self%on = .false.
597  else if (llc(i) <= self%rt_llc(i)) then
598  call self%boundaries%set (lbit(i))
599  self%mesh(i)%lower_boundary = .true.
600  end if
601  !
602  if ((llc(i) > self%rt_urc(i))) then
603  call self%patch%set(bits%no_rt)
604  self%on = .false.
605  else if (urc(i) >= self%rt_urc(i)) then
606  call self%boundaries%set (ubit(i))
607  self%mesh(i)%upper_boundary = .true.
608  end if
609  end do
610 END SUBROUTINE setup_boundaries
611 
612 !===============================================================================
613 !> This is called only for MHD tasks, just before they update, with self = main,
614 !> so this is suitable place to add the resulting net radiative heating to the
615 !> heating per unit volume.
616 !===============================================================================
617 SUBROUTINE pre_update (self)
618  class(rt_t):: self
619  !.............................................................................
620  class(rt_t), pointer:: omega
621  integer:: l(3), u(3), i_omega, i_bin
622  integer, save:: itimer=0
623  !-----------------------------------------------------------------------------
624  if (.not. self%on) &
625  return
626  call trace%begin ('rt_t%pre_update', itimer=itimer)
627  !---------------------------------------------------------------------------
628  ! -- Sum Q over angles and bin
629  self%qr = 0.0
630  ! -- multiply with the stored net intensities
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)
636  end do
637  end do
638  call self%patch%aux%register ('qr', self%qr)
639  !---------------------------------------------------------------------------
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
644  !---------------------------------------------------------------------------
645  call trace%end (itimer)
646 END SUBROUTINE pre_update
647 
648 !===============================================================================
649 !> This is called from MHD tasks, after updates, with self = main RT task.
650 !> Sets the time step of the RT tasks equal to the MHD timestep, which was
651 !> recently computed, and collects u_max and modifies the mhd%u_max
652 !===============================================================================
653 SUBROUTINE post_update (self)
654  class(rt_t) :: self
655  !-----------------------------------------------------------------------------
656  class(gpatch_t), pointer:: mhd
657  class(rt_t), pointer:: omega
658  integer:: i
659  !-----------------------------------------------------------------------------
660  mhd => self%patch
661  if (verbose == -1) then
662  !$omp atomic
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
667  !$omp atomic write
668  timer%n_mhd = 0
669  !$omp atomic write
670  timer%n_solve = 0
671  end if
672  end if
673  if (.not.self%on) return
674  call trace%begin ('rt_t%post_update')
675  !-----------------------------------------------------------------------------
676  ! We now know the timestep size and the time, and we can set these for all the
677  ! RT solvers, in such a way that they update in the desired order.
678  !-----------------------------------------------------------------------------
679  self%dtime = mhd%dtime
680  self%time = mhd%time + self%grace*self%dtime
681  do i=1,self%n_omega
682  omega => self%omega(i)
683  omega%time = mhd%time
684  !-------------------------------------------------------------------------
685  ! Inclined tasks may advance their time a factor mu slower, and thus forced
686  ! to update more often, to propagate the solution by about the same amount
687  ! vertically. NOTE (FIXME): Purely horizontal rays need special treatment
688  !-------------------------------------------------------------------------
689  if (self%vertical_propagation) then
690  omega%dtime = mhd%dtime*abs(omega%mu)
691  else
692  omega%dtime = mhd%dtime + 2*self%grace*self%dtime
693  end if
694  end do
695  !-----------------------------------------------------------------------------
696  call trace%end()
697 END SUBROUTINE post_update
698 
699 !===============================================================================
700 !> The rt_t%update procedure gets called directly by the task dispatcher, when
701 !> it has concluded that the task is ready to update. An update of the main RT
702 !> task consists of summing up the results from the RT sub-tasks, while the
703 !> sub-task updates perform the actual RT solutions along rays.
704 !===============================================================================
705 SUBROUTINE update (self)
706  class(rt_t):: self
707  !.............................................................................
708  class(rt_t), pointer :: omega
709  integer:: i_omega, i_bin
710  real, dimension(:,:,:), pointer:: q
711  real(8):: time
712  integer, save:: itimer=0
713  !-----------------------------------------------------------------------------
714  if (.not. self%on) &
715  return
716  call trace%begin ('rt_t%update', itimer=itimer)
717  !-----------------------------------------------------------------------------
718  ! The omega0 task is responsible for calling rt_t%rt_eos to compute the RT EOS
719  ! data. The nbor list has triggered download of nbor MHD data, before this call
720  !-----------------------------------------------------------------------------
721  if (self%i_omega==0) then
722  call self%rt_eos (self%time+self%dtime)
723  !-----------------------------------------------------------------------------
724  ! RT sub-task, solving for the radiation field in a given direction
725  !-----------------------------------------------------------------------------
726  else
727  ! -- Call RT solver, storing result in %new slot, at time = t(%it)
728  self%q => self%mem(:,:,:,:,self%new,1)
729  call self%solve ()
730  if (verbose == -1) then
731  !$omp atomic
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
735  end if
736  end if
737  call trace%end (itimer)
738 END SUBROUTINE update
739 
740 !===============================================================================
741 !> Compute the RT EOS quantities for the patch, using self%eos_time and an OMP
742 !> critical region to guard against duplicate calculations when multi-threading.
743 !===============================================================================
744 SUBROUTINE rt_eos (self, time)
745  USE units_mod
746  class(rt_t):: self
747  real(8):: time
748  !.............................................................................
749  class(gpatch_t), pointer :: patch
750  integer, save:: itimer=0
751  real, dimension(:,:,:), pointer:: d, e, pgas, fmax
752  !-----------------------------------------------------------------------------
753  if (time > self%eos_time) then
754  !$omp critical (eos_cr)
755  if (time > self%eos_time) then
756  call omp_eos
757  self%eos_time = time
758  if (verbose == -1) &
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
764  end if
765  !$omp end critical (eos_cr)
766  end if
767  return
768 contains
769 !-------------------------------------------------------------------------------
770 !-------------------------------------------------------------------------------
771 subroutine omp_eos
772  integer:: m(4), ib
773  real:: u_max, stefan
774  !-----------------------------------------------------------------------------
775  call trace%begin ('rt_t%rt_eos', itimer=itimer)
776  patch => self%patch
777  m = shape(self%rk)
778  allocate (pgas(m(1),m(2),m(3)),fmax(m(1),m(2),m(3)))
779  !-----------------------------------------------------------------------------
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')
785  end if
786  call eos%lookup_table (m, d=d, e=e, rk=self%rk, src=self%src, tt=self%tt, pg=pgas)
787  ! ----------------------------------------------------------------------------
788  ! Compute radiative diffusion speed, and set default omega%u_max
789  ! ----------------------------------------------------------------------------
790  self%u_max = 0.0
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))
794  !fmax = (pgas/d)**4/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)
798  end do
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')
805  end if
806  ! ----------------------------------------------------------------------------
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.
811  self%dt_fixed = -1.0
812  end if
813  end if
814  call trace%end (itimer)
815 end subroutine omp_eos
816 END SUBROUTINE rt_eos
817 
818 !===============================================================================
819 !> Call integral short characteristics RT solver.
820 !===============================================================================
821 SUBROUTINE solve (self)
822  class(rt_t):: self
823  !.............................................................................
824  integer, save:: itimer=0
825  !-----------------------------------------------------------------------------
826  if (self%i_omega==0) return
827  call trace%begin ('rt_t%solve', itimer=itimer)
828  !-----------------------------------------------------------------------------
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, &
832  self%rt_hat)
833  call self%patch%aux%register ('rkap', self%rk )
834  call self%patch%aux%register ('src', self%src)
835  call trace%end (itimer)
836 END SUBROUTINE solve
837 
838 !===============================================================================
839 !> For diagnostic output of tasks and nbor lists to data/run/rank_RRRRR.log
840 !===============================================================================
841 SUBROUTINE task_info (self)
842  class(rt_t):: self
843  !-----------------------------------------------------------------------------
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
850 
851 !===============================================================================
852 !> For diagnostic output of tasks and nbor lists to data/run/rank_RRRRR.log
853 !===============================================================================
854 SUBROUTINE nbor_info (self, nbor)
855  class(rt_t):: self
856  class(link_t):: nbor
857  !-----------------------------------------------------------------------------
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
864 
865 !===============================================================================
866 !> Intercept of courant_condition() call
867 !===============================================================================
868 SUBROUTINE courant_condition (self, detailed_timer)
869  class(rt_t):: self
870  logical, optional:: detailed_timer
871  !-----------------------------------------------------------------------------
872  call trace%begin ('rt_t%courant_condition')
873  self%patch%u_max = max(self%patch%u_max, self%u_max)
874  call trace%end()
875 END SUBROUTINE courant_condition
876 
877 !===============================================================================
878 !> Auxiliary output
879 !===============================================================================
880 SUBROUTINE output (self)
881  class(rt_t):: self
882  !-----------------------------------------------------------------------------
883 END SUBROUTINE output
884 
885 END MODULE rt_mod
Each thread uses a private timer data type, with arrays for start time and total time for each regist...
Definition: timer_mod.f90:11
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.
Definition: gpatch_mod.f90:4
Support tic/toc timing, as in MATLAB, and accurate wallclock() function. The timing is generally much...
Module with list handling for generic class task_t objects.
Definition: list_mod.f90:4
Define code units, in terms of (here) CGS units.
Definition: scaling_mod.f90:4
Module to set up angles and bins for an RT patch with only absorption. This means that the only extra...
Definition: rt_mod.f90:46
Equation of state module for any sort of tables, provided by a reader.
Definition: eos_mod.f90:4
Fundamental constants in CGS and SI units.
Definition: units_mod.f90:4
Boundary conditions for centered variables, which have the physical boundary in-between mesh points l...
Fundamental constants in CGS and SI units.
Definition: math_mod.f90:4
Template module for mesh.
Definition: mesh_mod.f90:4
This index file has slot indices for all solver, all initially equal to zero It is the responsibility...
Definition: index_mod.f90:7
Definition: io_mod.f90:4
Module with which one can register any number of pointers to real or integer arrays, and then output the contents of the links later.
Definition: aux_mod.f90:14
Template module for tasks.
Definition: task_mod.f90:4