20 procedure(solve1),
pointer:: solver => null()
24 procedure:: prepare_layers
26 procedure:: calculate_offsets
29 procedure,
nopass:: diagnostics
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
40 SUBROUTINE init (self)
43 logical,
save:: first_time=.true.
44 namelist /rt_integral_params/ verbose, solver, order, dtau_min, dtau_dif, &
45 dtau_max, detailed_timer
52 read (io_unit%input, rt_integral_params, iostat=iostat)
53 write (io_unit%output, rt_integral_params)
62 write (stderr,*) solver
64 call io%abort (
'rt_integral%init: wrong solver')
72 SUBROUTINE solve (self, k, s, q, mesh, dir, mu)
74 real,
dimension(:,:,:,:),
pointer:: k, s, q
75 class(
mesh_t),
pointer:: mesh(:)
76 real:: mu3, phi3, mu(3)
78 real,
dimension(:,:,:,:),
allocatable:: k3, s3, q3
81 integer:: m(4), loc(1), dir, i, i1, i2, i3, ib, l(3), u(3), ii(3)
82 integer,
save:: itimer=0
84 call trc_begin (
'rt_integral_t%solve', itimer=itimer)
88 dr = 0.5*mesh(dir)%d/abs(mu(dir))
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))
108 if (mu(dir) > 0.0)
then 109 i1=l(3); i2=u(3); i3=+1
111 i1=u(3); i2=l(3); i3=-1
114 write(io_unit%log,*)
'rt_integral_t%solve: mu, phi, dir, dr =', mu, dir, 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)
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)
127 call self%solver (k3, s3, q3, q(:,:,i,:), dr, l, u, debug)
130 deallocate (k3, s3, q3)
131 if (verbose > 1)
then 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)
140 call trc_end (itimer)
148 SUBROUTINE prepare_layers (self, mu, dir, mesh, in, out, i, i3)
152 class(
mesh_t),
pointer:: mesh(:)
153 real,
dimension(:,:,:,:):: in, out
154 integer:: o(2,2), l(2), u(2)
156 integer,
save:: itimer=0
158 call trc_begin (
'rt_integral_t%prepare_layers', itimer=itimer)
160 call self%calculate_offsets (mu, dir, mesh, o, p)
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)
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
188 SUBROUTINE calculate_offsets (self, mu, dir, mesh, o, p)
190 class(
mesh_t),
pointer:: mesh(:)
191 integer:: dir, o(2,2)
194 real(8),
pointer,
dimension(:):: x, y, z
197 integer,
save:: itimer=0
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 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
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
215 p(:,j) = (2*j-3)*dpdi
218 o = max(-1,min(0,floor(p)))
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 =', &
225 END SUBROUTINE calculate_offsets
230 SUBROUTINE bilinear (self, in, out, o, p, l, u)
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
237 integer:: i1, i2, ib, m(3)
238 integer,
save:: itimer=0
240 call trc_begin (
'rt_integral_t%bilinear', itimer=itimer)
246 if (all(abs(o+p) < 1e-5))
then 249 write(io_unit%log,*)
'rt_integral_t%bilinear: no interpolation' 253 else if (abs(o(1)+p(1)) < 1e-5)
then 255 write(stderr,*)
'rt_integral_t%bilinear: 2nd index, o, p =', &
257 call check_index_2 (
'(1)')
262 q(2)*in(i1,i2+o(2),ib) + p(2)*in(i1,i2+o(2)+1,ib)
269 else if (abs(o(2)+p(2)) < 1e-5)
then 271 write(stderr,*)
'rt_integral_t%bilinear: 1st index, o, p =', &
273 call check_index_1 (
'(1)')
278 q(1)*in(i1+o(1),i2,ib) + p(1)*in(i1+o(1)+1,i2,ib)
287 write(stderr,*)
'rt_integral_t%bilinear: o, p =', &
289 call check_index_1 (
'(2)')
290 call check_index_2 (
'(2)')
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)
303 call trc_end (itimer)
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')
314 end subroutine check_index_1
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')
324 end subroutine check_index_2
325 END SUBROUTINE bilinear
335 SUBROUTINE solve1 (self, k, s, qin, qout, dr, l, u, debug)
337 real,
dimension(:,:,:,:) :: k, s, qin
338 real,
dimension(:,:,:) :: qout
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
349 call trc_begin (
'rt_integral_t%solve1', itimer=itimer)
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)))
360 if (dtaumin > dtau_max)
then 368 else if (dtaumin > dtau_dif)
then 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
379 else if (dtaumax < dtau_min)
then 382 qout(i1,i2,ib) = qin(i1,i2,ib,1)*ex0
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 &
401 n_cases(icase) = n_cases(icase)+1
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)
413 call trc_end (itimer)
414 END SUBROUTINE solve1
416 SUBROUTINE solve2 (self, k, s, qin, qout, dr, l, u, debug)
418 real,
dimension(:,:,:,:) :: k, s, qin
419 real,
dimension(:,:,:) :: qout
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
430 call trc_begin (
'rt_integral_t%solve2', itimer=itimer)
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))
437 ex0 = exp(-dtau1(i1))
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))
445 qout(i1,i2,ib) = qin(i1,i2,ib,1)*ex0 &
447 + merge(d2sdtau2 *ex2, 0.0, dtau1(i1) > dtau_min)
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)
460 call trc_end (itimer)
461 END SUBROUTINE solve2
464 SUBROUTINE diagnostics (label)
465 character(len=*),
optional:: label
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)
472 write(io_unit%log,*)
'fractions of thin, full, diff, and zero cases: ', &
475 END SUBROUTINE diagnostics
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)
487 SUBROUTINE trc_end (itimer)
488 integer,
optional:: itimer
489 if (detailed_timer)
then 490 call timer%end (itimer)
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...
Template module for mesh.