DISPATCH
rt_integral_mod.f90
1 !===============================================================================
2 !> Given a direction Omega, which may align mainly with the x-. y-, or z-axis,
3 !> solve the integral form of the RT on short characteristics, extending across
4 !> only three layers perpendicular to the main direction.
5 !>
6 !> This is done by transposing the incoming problem, defined by rk, src, and q,
7 !> in connection with doing bilinear interpolation in the planes perpendicular
8 !> to the main RT axis.
9 !===============================================================================
11  USE timer_mod
12  USE io_mod
13  USE io_unit_mod
14  USE trace_mod
15  USE mesh_mod
16  USE omp_mod
17  implicit none
18  private
19  type, public:: rt_integral_t
20  procedure(solve1), pointer:: solver => null()
21  contains
22  procedure:: init
23  procedure:: solve
24  procedure:: prepare_layers
25  procedure:: bilinear
26  procedure:: calculate_offsets
27  procedure:: solve1
28  procedure:: solve2
29  procedure, nopass:: diagnostics
30  end type
31  integer, save:: verbose=0, order=3, solver=2
32  logical:: detailed_timer=.false.
33  real:: dtau_min=1e-3, dtau_dif=1e3, dtau_max=1e4
34  integer(8), save:: n_cases(4)=0
35  type(rt_integral_t), public:: rt_integral
36 CONTAINS
37 
38 !===============================================================================
39 !===============================================================================
40 SUBROUTINE init (self)
41  class(rt_integral_t):: self
42  integer:: iostat
43  logical, save:: first_time=.true.
44  namelist /rt_integral_params/ verbose, solver, order, dtau_min, dtau_dif, &
45  dtau_max, detailed_timer
46  !-----------------------------------------------------------------------------
47  if (first_time) then
48  !$omp critical (input_cr)
49  if (first_time) then
50  first_time=.false.
51  rewind(io_unit%input)
52  read (io_unit%input, rt_integral_params, iostat=iostat)
53  write (io_unit%output, rt_integral_params)
54  end if
55  !$omp end critical (input_cr)
56  select case (solver)
57  case (1)
58  self%solver => solve1
59  case (2)
60  self%solver => solve2
61  case default
62  write (stderr,*) solver
63  flush (stderr)
64  call io%abort ('rt_integral%init: wrong solver')
65  end select
66  end if
67 END SUBROUTINE init
68 
69 !===============================================================================
70 !> Solve the radiative transfer for all bins, along the closest axis direction
71 !===============================================================================
72 SUBROUTINE solve (self, k, s, q, mesh, dir, mu)
73  class(rt_integral_t):: self
74  real, dimension(:,:,:,:), pointer:: k, s, q
75  class(mesh_t), pointer:: mesh(:)
76  real:: mu3, phi3, mu(3)
77  !.............................................................................
78  real, dimension(:,:,:,:), allocatable:: k3, s3, q3
79  real:: dr
80  logical:: debug
81  integer:: m(4), loc(1), dir, i, i1, i2, i3, ib, l(3), u(3), ii(3)
82  integer, save:: itimer=0
83  !-----------------------------------------------------------------------------
84  call trc_begin ('rt_integral_t%solve', itimer=itimer)
85  debug = .false.
86  m = shape(k)
87  !-----------------------------------------------------------------------------
88  dr = 0.5*mesh(dir)%d/abs(mu(dir))
89  !-----------------------------------------------------------------------------
90  ! Transpose dimensions to solver index coordinates and allocate
91  !-----------------------------------------------------------------------------
92  if (dir==1) then
93  ii = [2,3,1]
94  else if (dir==2) then
95  ii = [3,1,2]
96  else
97  ii = [1,2,3]
98  end if
99  m(1:3) = m(ii)
100  allocate (k3(m(1),m(2),m(4),3))
101  allocate (s3(m(1),m(2),m(4),3))
102  allocate (q3(m(1),m(2),m(4),3))
103  !-----------------------------------------------------------------------------
104  ! Determine loop limits and loop direction
105  !-----------------------------------------------------------------------------
106  l = mesh(ii)%li
107  u = mesh(ii)%ui
108  if (mu(dir) > 0.0) then
109  i1=l(3); i2=u(3); i3=+1
110  else
111  i1=u(3); i2=l(3); i3=-1
112  end if
113  if (verbose > 1) &
114  write(io_unit%log,*) 'rt_integral_t%solve: mu, phi, dir, dr =', mu, dir, i1, i2, i3
115  !-----------------------------------------------------------------------------
116  ! Prepare layers for the single layer solver and call it
117  !-----------------------------------------------------------------------------
118  do i=i1,i2,i3
119  call self%prepare_layers (mu, dir, mesh, k, k3, i, i3)
120  call self%prepare_layers (mu, dir, mesh, s, s3, i, i3)
121  call self%prepare_layers (mu, dir, mesh, q, q3, i, i3)
122  if (dir==1) then
123  call self%solver (k3, s3, q3, q(i,:,:,:), dr, l, u, debug)
124  else if (dir==2) then
125  call self%solver (k3, s3, q3, q(:,i,:,:), dr, l, u, debug)
126  else
127  call self%solver (k3, s3, q3, q(:,:,i,:), dr, l, u, debug)
128  end if
129  end do
130  deallocate (k3, s3, q3)
131  if (verbose > 1) then
132  do ib=1,m(4)
133  write(io_unit%log,'(a,i4,1p,10e10.2)') ' rt_integral_t%solve: q =', &
134  ib, q(l(1):l(1)+9,(l(2)+u(2))/2,(l(3)+u(3))/2,ib)
135  end do
136  end if
137  if (verbose > 0) &
138  call diagnostics
139  !-----------------------------------------------------------------------------
140  call trc_end (itimer)
141 END SUBROUTINE solve
142 
143 !===============================================================================
144 !> Prepare three layers; one centered on the mesh points, one before and one
145 !> after, in the RT direction. The indices in the layers are (i1,i2,ib,1:3),
146 !> where i1-i3 are 3D coords, ib is bin index, and the last is layer index
147 !===============================================================================
148 SUBROUTINE prepare_layers (self, mu, dir, mesh, in, out, i, i3)
149  class(rt_integral_t):: self
150  real:: mu(3)
151  integer:: dir, i, i3
152  class(mesh_t), pointer:: mesh(:)
153  real, dimension(:,:,:,:):: in, out
154  integer:: o(2,2), l(2), u(2)
155  real:: p(2,2)
156  integer, save:: itimer=0
157  !-----------------------------------------------------------------------------
158  call trc_begin ('rt_integral_t%prepare_layers', itimer=itimer)
159  !-----------------------------------------------------------------------------
160  call self%calculate_offsets (mu, dir, mesh, o, p)
161  if (dir==1) then
162  l = [mesh(2)%li, mesh(3)%li]
163  u = [mesh(2)%ui, mesh(3)%ui]
164  call self%bilinear (in(i-i3,:,:,:), out(:,:,:,1), o(:,1), p(:,1), l, u)
165  out(:,:,:,2) = in(i,:,:,:)
166  call self%bilinear (in(i+i3,:,:,:), out(:,:,:,3), o(:,2), p(:,2), l, u)
167  else if (dir==2) then
168  l = [mesh(1)%li, mesh(3)%li]
169  u = [mesh(1)%ui, mesh(3)%ui]
170  call self%bilinear (in(:,i-i3,:,:), out(:,:,:,1), o(:,1), p(:,1), l, u)
171  out(:,:,:,2) = in(:,i,:,:)
172  call self%bilinear (in(:,i+i3,:,:), out(:,:,:,3), o(:,2), p(:,2), l, u)
173  else if (dir==3) then
174  l = [mesh(1)%li, mesh(2)%li]
175  u = [mesh(1)%ui, mesh(2)%ui]
176  call self%bilinear (in(:,:,i-i3,:), out(:,:,:,1), o(:,1), p(:,1), l, u)
177  out(:,:,:,2) = in(:,:,i,:)
178  call self%bilinear (in(:,:,i+i3,:), out(:,:,:,3), o(:,2), p(:,2), l, u)
179  end if
180  if (verbose > 1) &
181  write(io_unit%log,*) 'rt_integral_t%prepare_layers: dir, l, u =', dir, l, u
182  call trc_end (itimer)
183 END SUBROUTINE prepare_layers
184 
185 !===============================================================================
186 !> Calculate perpendicular index offsets from one layer to the next
187 !===============================================================================
188 SUBROUTINE calculate_offsets (self, mu, dir, mesh, o, p)
189  class(rt_integral_t):: self
190  class(mesh_t), pointer:: mesh(:)
191  integer:: dir, o(2,2)
192  real:: mu(3), p(2,2)
193  !.............................................................................
194  real(8), pointer, dimension(:):: x, y, z
195  real:: dpdi(2)
196  integer:: j
197  integer, save:: itimer=0
198  !-----------------------------------------------------------------------------
199  ! -- Fractional steps in perpendicular indices per parallel index
200  if (dir==1) then
201  ! -- when dir==1, the 1st index is y, the 2nd is z
202  dpdi(1) = mu(2)/abs(mu(1))*mesh(1)%d/mesh(2)%d
203  dpdi(2) = mu(3)/abs(mu(1))*mesh(1)%d/mesh(3)%d
204  else if (dir==2) then
205  ! -- when dir==2, the 1st index is z, the 2nd is x
206  dpdi(1) = mu(3)/abs(mu(2))*mesh(2)%d/mesh(3)%d
207  dpdi(2) = mu(1)/abs(mu(2))*mesh(2)%d/mesh(1)%d
208  else
209  ! -- when dir==3, the 1st index is x, the 2nd is y
210  dpdi(1) = mu(1)/abs(mu(3))*mesh(3)%d/mesh(1)%d
211  dpdi(2) = mu(2)/abs(mu(3))*mesh(3)%d/mesh(2)%d
212  end if
213  ! -- backward and forwards one index
214  do j=1,2
215  p(:,j) = (2*j-3)*dpdi
216  end do
217  ! -- prevent o=1 p=0.0 from happening
218  o = max(-1,min(0,floor(p)))
219  p = p - o
220  if (verbose > 1) then
221  write(io_unit%log,'(i2,2x,a,2x,3f7.3,2x,2f7.3,2(2x,2i3),0p,2(2x,2f7.3))') &
222  omp%thread, 'rt_integral_t%calculate_offsets: mu, dpdi, o, p =', &
223  mu, dpdi, o, p
224  end if
225 END SUBROUTINE calculate_offsets
226 
227 !===============================================================================
228 !> Bilinear interpolation in a plane
229 !===============================================================================
230 SUBROUTINE bilinear (self, in, out, o, p, l, u)
231  class(rt_integral_t):: self
232  real, dimension(:,:,:), intent(in) :: in
233  real, dimension(:,:,:), intent(out) :: out
234  real, dimension(2):: p, q
235  integer, dimension(2):: o, l, u
236  !.............................................................................
237  integer:: i1, i2, ib, m(3)
238  integer, save:: itimer=0
239  !-----------------------------------------------------------------------------
240  call trc_begin ('rt_integral_t%bilinear', itimer=itimer)
241  m = shape(in)
242  q = 1.0-p
243  !-----------------------------------------------------------------------------
244  ! Purely vertical case
245  !-----------------------------------------------------------------------------
246  if (all(abs(o+p) < 1e-5)) then
247  out = in
248  if (verbose > 1) &
249  write(io_unit%log,*) 'rt_integral_t%bilinear: no interpolation'
250  !-----------------------------------------------------------------------------
251  ! Interpolation only in y
252  !-----------------------------------------------------------------------------
253  else if (abs(o(1)+p(1)) < 1e-5) then
254  if (verbose > 1) &
255  write(stderr,*) 'rt_integral_t%bilinear: 2nd index, o, p =', &
256  o(2), p(2)
257  call check_index_2 ('(1)')
258  do ib=1,m(3)
259  do i2=l(2),u(2)
260  do i1=l(1),u(1)
261  out(i1,i2,ib) = &
262  q(2)*in(i1,i2+o(2),ib) + p(2)*in(i1,i2+o(2)+1,ib)
263  end do
264  end do
265  end do
266  !-----------------------------------------------------------------------------
267  ! Interpolation only in x
268  !-----------------------------------------------------------------------------
269  else if (abs(o(2)+p(2)) < 1e-5) then
270  if (verbose > 1) &
271  write(stderr,*) 'rt_integral_t%bilinear: 1st index, o, p =', &
272  o(1), p(1)
273  call check_index_1 ('(1)')
274  do ib=1,m(3)
275  do i2=l(2),u(2)
276  do i1=l(1),u(1)
277  out(i1,i2,ib) = &
278  q(1)*in(i1+o(1),i2,ib) + p(1)*in(i1+o(1)+1,i2,ib)
279  end do
280  end do
281  end do
282  !-----------------------------------------------------------------------------
283  ! Interpolation in x and y
284  !-----------------------------------------------------------------------------
285  else
286  if (verbose > 1) &
287  write(stderr,*) 'rt_integral_t%bilinear: o, p =', &
288  o, p
289  call check_index_1 ('(2)')
290  call check_index_2 ('(2)')
291  do ib=1,m(3)
292  do i2=l(2),u(2)
293  do i1=l(1),u(1)
294  out(i1,i2,ib) = &
295  q(1)*q(2)*in(i1+o(1) ,i2+o(2) ,ib) + &
296  p(1)*q(2)*in(i1+o(1)+1,i2+o(2) ,ib) + &
297  q(1)*p(2)*in(i1+o(1) ,i2+o(2)+1,ib) + &
298  p(1)*p(2)*in(i1+o(1)+1,i2+o(2)+1,ib)
299  end do
300  end do
301  end do
302  end if
303  call trc_end (itimer)
304 contains
305 !===============================================================================
306 subroutine check_index_1 (label)
307  character(len=*):: label
308  do i1=l(1),u(1),u(1)-l(1)
309  if (i1+o(1) < 1 .or. i1+o(1)+1 > size(in,1)) then
310  write(stderr,*) 'index 1 outside range '//label, i1, size(in,1), o(1), p(1)
311  call io%abort ('rt_integral_t%bilinear: index outside range')
312  end if
313  end do
314 end subroutine check_index_1
315 !===============================================================================
316 subroutine check_index_2 (label)
317  character(len=*):: label
318  do i2=l(2),u(2),u(2)-l(2)
319  if (i2+o(2) < 1 .or. i2+o(2)+1 > size(in,2)) then
320  write(stderr,*) 'index 2 outside range '//label, i2, size(in,2), o(2), p(2)
321  call io%abort ('rt_integral_t%bilinear: index outside range')
322  end if
323  end do
324 end subroutine check_index_2
325 END SUBROUTINE bilinear
326 
327 !===============================================================================
328 !> Integral solver for one layer. The exact integral solution Q=I-S for a 2nd
329 !> order polynomial has coefficient for the upstream Q, and for the 1st and 2nd
330 !> derivative of S. However, if one is doing a solution both forwards and
331 !> backward, the term in the 1st derivative cancels.
332 !>
333 !> The order of indices is in general permuted in the call, with
334 !===============================================================================
335 SUBROUTINE solve1 (self, k, s, qin, qout, dr, l, u, debug)
336  class(rt_integral_t):: self
337  real, dimension(:,:,:,:) :: k, s, qin
338  real, dimension(:,:,:) :: qout
339  integer:: l(3), u(3)
340  real:: dr
341  logical:: debug
342  !.............................................................................
343  integer:: i1, i2, ib, m(4), icase
344  real:: dsdtau1, dsdtau2, d2sdtau2, ex0, ex2, dtaumin, dtaumax
345  real(8), pointer:: z(:)
346  real, dimension(size(k,1)):: dtau1, dtau2
347  integer, save:: itimer=0
348  !-----------------------------------------------------------------------------
349  call trc_begin ('rt_integral_t%solve1', itimer=itimer)
350  m = shape(k)
351  do ib=1,m(3)
352  do i2=l(2),u(2)
353  dtau1(l(1):u(1)) = dr*(k(l(1):u(1),i2,ib,2)+k(l(1):u(1),i2,ib,1))
354  dtau2(l(1):u(1)) = dr*(k(l(1):u(1),i2,ib,2)+k(l(1):u(1),i2,ib,3))
355  dtaumin = minval(dtau1(l(1):u(1)))
356  dtaumax = maxval(dtau1(l(1):u(1)))
357  !-------------------------------------------------------------------------
358  ! For optical depth increments larger than dtau_max for all points: Q=0.0
359  !-------------------------------------------------------------------------
360  if (dtaumin > dtau_max) then
361  do i1=l(1),u(1)
362  qout(i1,i2,ib) = 0.0
363  end do
364  icase = 4
365  !-------------------------------------------------------------------------
366  ! For optical depth increments larger than dtau_dif: Q = diffusion approx
367  !-------------------------------------------------------------------------
368  else if (dtaumin > dtau_dif) then
369  do i1=l(1),u(1)
370  dsdtau1 = (s(i1,i2,ib,2) - s(i1,i2,ib,1))/dtau1(i1)
371  dsdtau2 = (s(i1,i2,ib,3) - s(i1,i2,ib,2))/dtau2(i1)
372  d2sdtau2 = (dsdtau2 - dsdtau1)*2.0/(dtau1(i1)+dtau2(i1))
373  qout(i1,i2,ib) = d2sdtau2
374  end do
375  icase = 3
376  !-------------------------------------------------------------------------
377  ! For optical depth increments all maller than dtau_min: linear absorption
378  !-------------------------------------------------------------------------
379  else if (dtaumax < dtau_min) then
380  do i1=l(1),u(1)
381  ex0 = 1.0-dtau1(i1)
382  qout(i1,i2,ib) = qin(i1,i2,ib,1)*ex0
383  end do
384  icase = 1
385  !-------------------------------------------------------------------------
386  ! For all other cases, use the full expression (assume sum of +/- directions)
387  !-------------------------------------------------------------------------
388  else
389  do i1=l(1),u(1)
390  dsdtau1 = (s(i1,i2,ib,2) - s(i1,i2,ib,1))/dtau1(i1)
391  dsdtau2 = (s(i1,i2,ib,3) - s(i1,i2,ib,2))/dtau2(i1)
392  d2sdtau2 = (dsdtau2 - dsdtau1)*2.0/(dtau1(i1)+dtau2(i1))
393  ex0 = exp(-dtau1(i1))
394  ex2 = (1.0-ex0) - dtau1(i1)*ex0
395  qout(i1,i2,ib) = qin(i1,i2,ib,1)*ex0 &
396  + d2sdtau2 *ex2
397  end do
398  icase = 2
399  end if
400  !$omp atomic
401  n_cases(icase) = n_cases(icase)+1
402  end do
403  if (verbose > 1 .and. debug) then
404  write(io_unit%log,1) ' rt_integral_t%solve1: ib, dt =', &
405  ib, icase, dtau1(l(1):l(1)+9)
406  write(io_unit%log,1) ' rt_integral_t%solve1: ib, qi =', &
407  ib, icase, qin(l(1):l(1)+9,(l(2)+u(2))/2,ib,1)
408  write(io_unit%log,1) ' rt_integral_t%solve1: ib, qo =', &
409  ib, icase, qout(l(1):l(1)+9,(l(2)+u(2))/2,ib)
410  1 format(a,i4,i2,1p,10e10.2)
411  end if
412  end do
413  call trc_end (itimer)
414 END SUBROUTINE solve1
415 
416 SUBROUTINE solve2 (self, k, s, qin, qout, dr, l, u, debug)
417  class(rt_integral_t):: self
418  real, dimension(:,:,:,:) :: k, s, qin
419  real, dimension(:,:,:) :: qout
420  integer:: l(3), u(3)
421  real:: dr
422  logical:: debug
423  !.............................................................................
424  integer:: i1, i2, ib, m(4), icase
425  real:: dsdtau1, dsdtau2, d2sdtau1, d2sdtau2, ex0, ex1, ex2, dtaumin, dtaumax
426  real(8), pointer:: z(:)
427  real, dimension(size(k,1)):: dtau1, dtau2
428  integer, save:: itimer=0
429  !-----------------------------------------------------------------------------
430  call trc_begin ('rt_integral_t%solve2', itimer=itimer)
431  m = shape(k)
432  do ib=1,m(3)
433  do i2=l(2),u(2)
434  dtau1(l(1):u(1)) = dr*(k(l(1):u(1),i2,ib,2)+k(l(1):u(1),i2,ib,1))
435  dtau2(l(1):u(1)) = dr*(k(l(1):u(1),i2,ib,2)+k(l(1):u(1),i2,ib,3))
436  do i1=l(1),u(1)
437  ex0 = exp(-dtau1(i1))
438  ex1 = 1.-ex0
439  ex2 = ex1 - dtau1(i1)*ex0
440  dsdtau1 = (s(i1,i2,ib,2) - s(i1,i2,ib,1))/dtau1(i1)
441  dsdtau2 = (s(i1,i2,ib,3) - s(i1,i2,ib,2))/dtau2(i1)
442  d2sdtau1 = (dsdtau2 * dtau1(i1) + dsdtau1 * dtau2(i1))/(dtau1(i1)+dtau2(i1))
443  d2sdtau2 = (dsdtau2 - dsdtau1)*2.0/(dtau1(i1)+dtau2(i1))
444 
445  qout(i1,i2,ib) = qin(i1,i2,ib,1)*ex0 &
446  - d2sdtau1 *ex1 &
447  + merge(d2sdtau2 *ex2, 0.0, dtau1(i1) > dtau_min)
448  end do
449  end do
450  if (verbose > 1 .and. debug) then
451  write(io_unit%log,1) ' rt_integral_t%solve2: ib, dt =', &
452  ib, icase, dtau1(l(1):l(1)+9)
453  write(io_unit%log,1) ' rt_integral_t%solve2: ib, qi =', &
454  ib, icase, qin(l(1):l(1)+9,(l(2)+u(2))/2,ib,1)
455  write(io_unit%log,1) ' rt_integral_t%solve2: ib, qo =', &
456  ib, icase, qout(l(1):l(1)+9,(l(2)+u(2))/2,ib)
457  1 format(a,i4,i2,1p,10e10.2)
458  end if
459  end do
460  call trc_end (itimer)
461 END SUBROUTINE solve2
462 
463 !===============================================================================
464 SUBROUTINE diagnostics (label)
465  character(len=*), optional:: label
466  real:: f
467  f = 1./max(sum(n_cases),1_8)
468  if (present(label)) then
469  write(io_unit%log,*) 'fractions of thin, full, diff, and zero cases: ', &
470  n_cases*f, trim(label)
471  else
472  write(io_unit%log,*) 'fractions of thin, full, diff, and zero cases: ', &
473  n_cases*f
474  end if
475 END SUBROUTINE diagnostics
476 
477 !===============================================================================
478 SUBROUTINE trc_begin (label, itimer)
479  character(len=*):: label
480  integer, optional:: itimer
481  if (detailed_timer) then
482  call timer%begin (label, itimer)
483  end if
484 END SUBROUTINE
485 
486 !===============================================================================
487 SUBROUTINE trc_end (itimer)
488  integer, optional:: itimer
489  if (detailed_timer) then
490  call timer%end (itimer)
491  end if
492 END SUBROUTINE
493 
494 END MODULE rt_integral_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...
Template module for mesh.
Definition: mesh_mod.f90:4
Definition: io_mod.f90:4