9 #include "sll_memory.h"
10 #include "sll_working_precision.h"
28 sll_int32,
dimension(:, :, :),
pointer :: index
50 sll_real64,
dimension(:),
intent(in) ::
data
51 sll_real64,
dimension(:),
intent(in) :: eta
53 if (interpolator%boundary == 0)
then
67 sll_real64,
dimension(:),
intent(inout) :: data_in
68 sll_real64,
dimension(:),
intent(out) :: data_out
69 sll_int32,
dimension(:),
intent(in) :: dorder
70 sll_real64,
intent(in) ::displacement
71 logical,
intent(in) :: hiera
73 sll_int32 :: i1, i2, i3, k2, k3, counter, j
74 sll_int32,
dimension(3) :: l, no, ind_order
75 sll_int32,
dimension(:, :),
allocatable :: ind
77 sll_allocate(ind(interpolator%max_level + 1, 3), i1);
78 ind_order(dorder(1)) = 0
80 do i3 = 0, interpolator%levels(dorder(3))
81 ind_order(dorder(3)) = i3
83 no(dorder(3)) = max(2**(i3 - 1), 1);
84 do i2 = 0, min(interpolator%max_level - i3, interpolator%levels(dorder(2)))
85 ind_order(dorder(2)) = i2
86 no(dorder(2)) = max(2**(i2 - 1), 1);
87 ind(1, dorder(1)) = 0;
88 do k2 = 0, no(dorder(2)) - 1
89 ind(ind_order(dorder(2)) + 1, dorder(2)) = k2;
90 do k3 = 0, no(dorder(3)) - 1
91 ind(ind_order(dorder(3)) + 1, dorder(3)) = k3;
92 counter = interpolator%index( &
93 ind_order(1), ind_order(2), ind_order(3)) + &
94 ind(ind_order(1) + 1, 1)*no(2)*no(3) &
95 + ind(ind_order(2) + 1, 2)*no(3) + &
96 ind(ind_order(3) + 1, 3)
98 call interpolator%interpolate_disp_1d_periodic &
99 (displacement, dorder(1), &
100 min(interpolator%levels(dorder(1)), &
101 interpolator%max_level - ind_order(dorder(2)) - &
102 ind_order(dorder(3))), counter, data_in, data_out, hiera)
108 if (hiera .EQV. .false.)
then
110 do j = interpolator%order, 2, -1
111 call interpolator%dehierarchical_part_order &
113 interpolator%dim, 2, dorder, j)
116 call interpolator%dehierarchical_part(data_out, &
117 interpolator%dim, 2, dorder)
126 sll_int32 :: j, l1, l2, l3, level
128 sll_real64,
dimension(:),
intent(in) ::
data, eta
129 sll_real64,
dimension(3) :: eta_norm
130 sll_real64,
dimension(3) :: phi
131 sll_int32,
dimension(3) :: no, l
132 sll_int32,
dimension(:, :),
allocatable :: ind
136 sll_allocate(ind(0:interpolator%max_level, 1:3), j)
141 do j = 1, interpolator%dim
142 eta_norm(j) = (eta(j) - interpolator%eta_min(j))/interpolator%length(j)
143 eta_norm(j) = modulo(eta_norm(j), 1.0_f64)
146 do level = 2, interpolator%max_level
147 ind(level, j) = ind(level - 1, j)*2
148 if (eta_norm(j) > scale*real(ind(level, j) + 1, f64))
then
149 ind(level, j) = ind(level, j) + 1
151 scale = scale*0.5_f64
155 do level = 0, interpolator%max_level
156 do l1 = 0, min(level, interpolator%levels(1))
158 no(1) = max(2**(l1 - 1), 1)
159 do l2 = max(0, level - l1 - interpolator%levels(3)), min(level - l1, interpolator%levels(2))
161 no(2) = max(2**(l2 - 1), 1)
162 l(3) = level - l1 - l2
164 no(3) = max(2**(l(3) - 1), 1)
166 index = interpolator%index(l1, l2, l3) + ind(l1, 1)*no(2)*no(3) &
167 + ind(l2, 2)*no(3) + ind(l3, 3)
169 call interpolator%basis_function(real(2**(max(l(j), 1)), f64)*eta_norm(j) &
170 - real(2*ind(l(j), j), f64) - 1.0_f64, phi(j), &
171 interpolator%hierarchy(index)%function_type(j))
173 val = val +
data(index)*phi(1)*phi(2)*phi(3)
183 sll_int32 :: j, l1, l2, l3, level
185 sll_real64,
dimension(:),
intent(in) ::
data
186 sll_real64,
dimension(:),
intent(in) :: eta
187 sll_real64,
dimension(3) :: eta_norm
188 sll_real64,
dimension(3) :: phi
189 sll_int32,
dimension(3) :: no, l
190 sll_int32,
dimension(:, :),
allocatable :: ind
191 sll_real64 :: scale, phi1a, phi2a, phi3a
194 sll_allocate(ind(0:interpolator%max_level, 1:interpolator%dim), j)
197 ind(0:1, 1:interpolator%dim) = 0
199 do j = 1, interpolator%dim
200 eta_norm(j) = (eta(j) - interpolator%eta_min(j))/interpolator%length(j)
203 do level = 2, interpolator%levels(j)
204 ind(level, j) = ind(level - 1, j)*2
205 if (eta_norm(j) > scale*real(ind(level, j) + 1, f64))
then
206 ind(level, j) = ind(level, j) + 1
208 scale = scale*0.5_f64
212 do level = 0, interpolator%max_level
213 do l1 = 0, min(level, interpolator%levels(1))
220 do l2 = max(0, level - l1 - interpolator%levels(3)), min(level - l1, interpolator%levels(2))
225 no(2) = max(2**(l2 - 1), 1);
227 l3 = level - l1 - l2;
231 no(3) = max(2**(l3 - 1), 1);
236 index = interpolator%index(l1, l2, l3) + &
237 ind(l1, 1)*no(2)*no(3) &
238 + ind(l2, 2)*no(3) + ind(l3, 3)
240 do j = 1, interpolator%dim
241 call interpolator%basis_function(real(2**l(j), f64)*eta_norm(j) &
242 - real(2*ind(l(j), j), f64) - 1.0_f64, phi(j), &
243 interpolator%hierarchy(index)%function_type(j))
245 val = val +
data(index) &
246 *phi(1)*phi(2)*phi(3)
248 index = interpolator%index(l1, l2, l3) + &
249 ind(l1, 1)*no(2)*no(3) &
250 + ind(l2, 2)*no(3) + ind(l3, 3)
251 call interpolator%basis_function(real(2**l(1), f64)*eta_norm(1) &
252 - real(2*ind(l(1), 1), f64) - 1.0_f64, phi(1), &
253 interpolator%hierarchy(index)%function_type(1))
254 call interpolator%basis_function(real(2**l(2), f64)*eta_norm(2) &
255 - real(2*ind(l(2), 2), f64) - 1.0_f64, phi(2), &
256 interpolator%hierarchy(index)%function_type(2))
257 call interpolator%basis_function(eta_norm(3), phi(3), -1)
258 val = val +
data(index)*phi(1)*phi(2)*phi(3);
259 call interpolator%basis_function(eta_norm(3) - 1.0_f64, phi(3), -1)
260 val = val +
data(index + 1)*phi(1)*phi(2)*phi(3);
264 index = interpolator%index(l1, l2, l3) + &
265 ind(l1, 1)*no(2)*no(3) &
266 + ind(l2, 2)*no(3) + ind(l3, 3)
267 call interpolator%basis_function(real(2**l(1), f64)*eta_norm(1) &
268 - real(2*ind(l(1), 1), f64) - 1.0_f64, phi(1), &
269 interpolator%hierarchy(index)%function_type(1))
270 call interpolator%basis_function(real(2**l(3), f64)*eta_norm(3) &
271 - real(2*ind(l(3), 3), f64) - 1.0_f64, phi(3), &
272 interpolator%hierarchy(index)%function_type(3))
273 call interpolator%basis_function(eta_norm(2), phi(2), -1)
274 val = val +
data(index)*phi(1)*phi(2)*phi(3);
275 call interpolator%basis_function(eta_norm(2) - 1.0_f64, phi(2), -1)
276 val = val +
data(index + no(3))*phi(1)*phi(2)*phi(3);
278 index = interpolator%index(l1, l2, l3) + &
279 ind(l1, 1)*no(2)*no(3) + ind(l2, 2)*no(3) + ind(l3, 3);
280 call interpolator%basis_function(real(2**l(1), f64)*eta_norm(1) &
281 - real(2*ind(l(1), 1), f64) - 1.0_f64, phi(1), &
282 interpolator%hierarchy(index)%function_type(1))
283 call interpolator%basis_function(eta_norm(3), phi(3), -1)
284 call interpolator%basis_function(eta_norm(3) - 1.0_f64, phi3a, -1)
285 call interpolator%basis_function(eta_norm(2), phi(2), -1)
286 call interpolator%basis_function(eta_norm(2) - 1.0_f64, phi2a, -1)
287 val = val +
data(index)*phi(1)*phi(2)*phi(3) + &
288 data(index + 1)*phi(1)*phi(2)*phi3a + &
289 data(index + 2)*phi(1)*phi2a*phi(3) + &
290 data(index + 3)*phi(1)*phi2a*phi3a
296 index = interpolator%index(l1, l2, l3) + &
297 ind(l1, 1)*no(2)*no(3) &
298 + ind(l2, 2)*no(3) + ind(l3, 3)
300 do j = 2, interpolator%dim
301 call interpolator%basis_function(real(2**l(j), f64)*eta_norm(j) &
302 - real(2*ind(l(j), j), f64) - 1.0_f64, phi(j), &
303 interpolator%hierarchy(index)%function_type(j))
305 call interpolator%basis_function(eta_norm(1), phi(1), -1)
306 val = val +
data(index)*phi(1)*phi(2)*phi(3);
307 call interpolator%basis_function(eta_norm(1) - 1.0_f64, phi(1), -1)
308 val = val +
data(index + no(2)*no(3))*phi(1)*phi(2)*phi(3)
310 index = interpolator%index(l1, l2, l3) + &
311 ind(l1, 1)*no(2)*no(3) &
312 + ind(l2, 2)*no(3) + ind(l3, 3)
313 call interpolator%basis_function(eta_norm(1), phi(1), -1)
314 call interpolator%basis_function(eta_norm(1) - 1.0_f64, phi1a, -1)
315 call interpolator%basis_function(real(2**l(2), f64)*eta_norm(2) &
316 - real(2*ind(l(2), 2), f64) - 1.0_f64, phi(2), &
317 interpolator%hierarchy(index)%function_type(2))
318 call interpolator%basis_function(eta_norm(3), phi(3), -1)
319 call interpolator%basis_function(eta_norm(3) - 1.0_f64, phi3a, -1)
320 val = val +
data(index)*phi(1)*phi(2)*phi(3) + &
321 data(index + 1)*phi(1)*phi(2)*phi3a + &
322 data(index + no(2)*no(3))*phi1a*phi(2)*phi(3) + &
323 data(index + no(2)*no(3) + 1)*phi1a*phi(2)*phi3a
327 index = interpolator%index(l1, l2, l3) + &
328 ind(l1, 1)*no(2)*no(3) &
329 + ind(l2, 2)*no(3) + ind(l3, 3)
330 call interpolator%basis_function(eta_norm(1), phi(1), -1)
331 call interpolator%basis_function(eta_norm(1) - 1.0_f64, phi1a, -1)
332 call interpolator%basis_function(real(2**l(3), f64)*eta_norm(3) &
333 - real(2*ind(l(3), 3), f64) - 1.0_f64, phi(3), &
334 interpolator%hierarchy(index)%function_type(3))
335 call interpolator%basis_function(eta_norm(2), phi(2), -1)
336 call interpolator%basis_function(eta_norm(2) - 1.0_f64, phi2a, -1)
337 val = val +
data(index)*phi(1)*phi(2)*phi(3) + &
338 data(index + no(3))*phi(1)*phi2a*phi(3) + &
339 data(index + no(2)*no(3))*phi1a*phi(2)*phi(3) + &
340 data(index + no(2)*no(3) + no(3))*phi1a*phi2a*phi(3);
342 index = interpolator%index(l1, l2, l3) + &
343 ind(l1, 1)*no(2)*no(3) + ind(l2, 2)*no(3) + ind(l3, 3);
344 call interpolator%basis_function(eta_norm(1), phi(1), -1)
345 call interpolator%basis_function(eta_norm(1) - 1.0_f64, phi1a, -1)
346 call interpolator%basis_function(eta_norm(3), phi(3), -1)
347 call interpolator%basis_function(eta_norm(3) - 1.0_f64, phi3a, -1)
348 call interpolator%basis_function(eta_norm(2), phi(2), -1)
349 call interpolator%basis_function(eta_norm(2) - 1.0_f64, phi2a, -1)
350 val = val +
data(index)*phi(1)*phi(2)*phi(3) + &
351 data(index + 1)*phi(1)*phi(2)*phi3a + &
352 data(index + 2)*phi(1)*phi2a*phi(3) + &
353 data(index + 3)*phi(1)*phi2a*phi3a + &
354 data(index + 4)*phi1a*phi(2)*phi(3) + &
355 data(index + 5)*phi1a*phi(2)*phi3a + &
356 data(index + 6)*phi1a*phi2a*phi(3) + &
357 data(index + 7)*phi1a*phi2a*phi3a
375 sll_int32,
intent(in) :: dim
376 sll_real64,
intent(in) :: displacment_in
377 sll_comp64,
dimension(:),
intent(inout) :: data_in
378 sll_real64,
dimension(:),
intent(out) :: data_out
379 sll_real64:: displacement
381 displacement = displacment_in*2.0_f64*
sll_p_pi/interpolator%length(dim)
382 call displace(interpolator, dim, displacement, data_in);
383 call ispfft(interpolator, data_in, data_out);
398 interpolation_type, &
405 sll_real64,
dimension(:),
intent(in) :: eta_min
406 sll_real64,
dimension(:),
intent(in) :: eta_max
407 sll_int32,
dimension(:),
intent(in) :: levels
409 sll_int32,
intent(in) :: order
410 sll_int32,
intent(in) :: interpolation
411 sll_int32,
intent(in) :: interpolation_type
412 sll_int32,
intent(in) :: modified
417 sll_int32,
intent(in) :: boundary
420 sll_int32 :: i, j, k1, k2, k3, l1, l2, l3, l, counter
421 sll_int32 :: ierr, no1, no2, no3
422 sll_int32,
dimension(:),
allocatable :: novec, lvec, kvec
424 interpolator%dim = 3;
425 interpolator%modified = modified;
426 interpolator%boundary = boundary;
427 sll_allocate(lvec(interpolator%dim), ierr);
428 sll_allocate(kvec(interpolator%dim), ierr);
429 sll_allocate(novec(interpolator%dim), ierr);
430 interpolator%max_level = levels(1);
431 do l = 2, interpolator%dim
432 interpolator%max_level = max(levels(l), interpolator%max_level);
434 if (interpolator%modified == 1)
then
435 interpolator%max_level = interpolator%max_level + 2;
438 interpolator%size_basis = 0;
439 if (interpolator%boundary == 0)
then
440 do l = 0, interpolator%max_level
441 do l1 = 0, min(l, levels(1))
442 do l2 = max(0, l - l1 - levels(3)), min(l - l1, levels(2))
444 interpolator%size_basis = &
445 interpolator%size_basis + &
446 max(2**(l1 - 1), 1)*max(2**(l2 - 1), 1)*max(2**(l3 - 1), 1)
451 do l = 0, interpolator%max_level
452 do l1 = 0, min(l, levels(1))
453 do l2 = max(0, l - l1 - levels(3)), min(l - l1, levels(2))
470 interpolator%size_basis = &
471 interpolator%size_basis + no1*no2*no3;
477 call interpolator%initialize_sg(levels, order, interpolation, &
478 interpolation_type, eta_min, eta_max);
479 sll_allocate(interpolator%index(0:interpolator%levels(1), 0:interpolator%levels(2), 0:interpolator%levels(3)), ierr)
483 do l = 0, interpolator%max_level
484 interpolator%level_mapping(l) = counter;
485 do l1 = 0, min(l, interpolator%levels(1))
486 novec(1) = max(2**(l1 - 1), 1)
487 if (interpolator%boundary == 1 .AND. l1 == 0)
then
491 do l2 = max(0, l - l1 - interpolator%levels(3)), min(l - l1, interpolator%levels(2))
492 novec(2) = max(2**(l2 - 1), 1)
493 if (interpolator%boundary == 1 .AND. l2 == 0)
then
497 lvec(3) = l - l1 - l2;
499 novec(3) = max(2**(l3 - 1), 1)
500 if (interpolator%boundary == 1 .AND. l3 == 0)
then
504 interpolator%index(l1, l2, l3) = counter
505 do k1 = 0, novec(1) - 1
507 do k2 = 0, novec(2) - 1
509 do k3 = 0, novec(3) - 1
511 do j = 1, interpolator%dim
512 if (interpolator%boundary == 0)
then
517 (interpolator, counter, j, lvec, kvec, novec);
520 counter = counter + 1;
527 interpolator%level_mapping(interpolator%max_level + 1) = counter;
529 do i = 1, interpolator%size_basis
530 do j = 1, interpolator%dim
531 interpolator%hierarchy(i)%coordinate(j) = &
532 interpolator%hierarchy(i)%coordinate(j)* &
533 interpolator%length(j) + &
534 interpolator%eta_min(j)
546 sll_int32,
intent(in) :: cdim
547 sll_int32,
intent(in) :: counter
548 sll_int32,
dimension(:),
intent(in) :: lvecin
549 sll_int32,
dimension(:),
intent(in) :: kvecin
550 sll_int32,
dimension(:),
intent(in) :: novecin
552 sll_int32,
dimension(3) :: lvec, kvec, novec
553 sll_int32 :: jj, stride
555 do jj = 1, interpolator%dim
556 lvec(jj) = lvecin(jj);
557 kvec(jj) = kvecin(jj);
558 novec(jj) = novecin(jj);
562 interpolator%hierarchy(counter)%level(cdim) = ld;
563 interpolator%hierarchy(counter)%index_on_level(cdim) = kd;
566 interpolator%hierarchy(counter)%coordinate(cdim) = 0.0_f64
567 interpolator%hierarchy(counter)%parent(stride) = &
569 interpolator%hierarchy(counter)%parent(stride + 1) = &
571 interpolator%hierarchy(counter)%function_type(cdim) = 0
573 interpolator%hierarchy(counter)%coordinate(cdim) = &
574 1.0_f64/(2.0_f64**ld) + real(kd, f64)*1.0_f64/(2.0_f64**(ld - 1))
576 lvec(cdim) = lvec(cdim) - 1;
577 novec(cdim) = max(novec(cdim)/2, 1);
579 kvec(cdim) = modulo((kd + 1)/2, max(2**(ld - 2), 1));
580 interpolator%hierarchy(counter)%parent( &
581 modulo(kd, 2) + stride) = &
582 interpolator%hierarchy( &
583 interpolator%index(lvec(1), lvec(2), lvec(3)) + &
584 kvec(1)*novec(2)*novec(3) + &
586 kvec(3))%parent(stride)
589 interpolator%hierarchy(counter)%parent( &
590 modulo(kd + 1, 2) + stride) = &
591 interpolator%index(lvec(1), lvec(2), lvec(3)) + &
592 kvec(1)*novec(2)*novec(3) + kvec(2)*novec(3) + &
595 interpolator%hierarchy( &
596 interpolator%hierarchy(counter)%parent( &
597 modulo(kd + 1, 2) + stride))%children(modulo(kd, 2) + stride) = counter
599 interpolator%hierarchy( &
600 interpolator%hierarchy(counter)%parent( &
601 modulo(kd + 1, 2) + stride))%children(modulo(kd + 1, 2) + stride) = counter
603 if (interpolator%order == 1)
then
604 interpolator%hierarchy(counter)% &
605 function_type(cdim) = -1
606 elseif (ld == 1 .OR. interpolator%order == 2)
then
607 interpolator%hierarchy(counter)% &
608 function_type(cdim) = 1 + modulo(kd, 2)
609 elseif (ld == 2 .OR. interpolator%order == 3)
then
610 interpolator%hierarchy(counter)% &
611 function_type(cdim) = 3 + modulo(kd, 4)
613 interpolator%hierarchy(counter)% &
614 function_type(cdim) = 7 + modulo(kd, 8)
625 sll_int32,
intent(in) :: cdim
626 sll_int32,
intent(in) :: counter
627 sll_int32,
dimension(:),
intent(in) :: lvecin
628 sll_int32,
dimension(:),
intent(in) :: kvecin
629 sll_int32,
dimension(:),
intent(in) :: novecin
631 sll_int32,
dimension(:),
allocatable :: lvec, kvec, novec
632 sll_int32 :: jj, stride
634 sll_allocate(lvec(interpolator%dim), jj);
635 sll_allocate(kvec(interpolator%dim), jj);
636 sll_allocate(novec(interpolator%dim), jj);
637 do jj = 1, interpolator%dim
638 lvec(jj) = lvecin(jj);
639 kvec(jj) = kvecin(jj);
640 novec(jj) = novecin(jj);
644 interpolator%hierarchy(counter)%level(cdim) = ld;
647 interpolator%hierarchy(counter)%coordinate(cdim) = real(kd, f64);
648 interpolator%hierarchy(counter)%parent(stride) = &
650 interpolator%hierarchy(counter)%parent(stride + 1) = &
652 interpolator%hierarchy(counter)%function_type(cdim) = -1
654 interpolator%hierarchy(counter)%coordinate(cdim) = &
655 1.0_f64/(2.0_f64**ld) + real(kd, f64)*1.0_f64/(2.0_f64**(ld - 1))
657 lvec(cdim) = lvec(cdim) - 1;
658 novec(cdim) = novec(cdim)/2;
662 interpolator%hierarchy(counter)%parent(stride) = &
663 interpolator%index(lvec(1), lvec(2), lvec(3)) + &
664 kvec(1)*novec(2)*novec(3) + kvec(2)*novec(3) + kvec(3)
666 interpolator%hierarchy(counter)%parent( &
667 stride + 1) = interpolator%hierarchy(counter)%parent( &
668 stride) + novec(2)*novec(3)
669 elseif (cdim == 2)
then
670 interpolator%hierarchy(counter)%parent( &
671 stride + 1) = interpolator%hierarchy(counter)%parent( &
674 interpolator%hierarchy(counter)%parent( &
675 stride + 1) = interpolator%hierarchy(counter)%parent( &
678 interpolator%hierarchy(interpolator%hierarchy(counter)%parent(stride))%children &
684 kvec(cdim) = modulo((kd + 1)/2, max(2**(ld - 2), 1));
685 kvec(cdim) = (kd + 1)/2;
686 if (kvec(cdim) < novec(cdim))
then
687 interpolator%hierarchy(counter)%parent( &
688 modulo(kd, 2) + stride) = &
689 interpolator%hierarchy( &
690 interpolator%index(lvec(1), lvec(2), lvec(3)) + &
691 kvec(1)*novec(2)*novec(3) + kvec(2)*novec(3) + kvec(3))% &
694 kvec(cdim) = kvec(cdim) - 1;
695 interpolator%hierarchy(counter)%parent( &
696 modulo(kd, 2) + stride) = &
697 interpolator%hierarchy( &
698 interpolator%index(lvec(1), lvec(2), lvec(3)) + &
699 kvec(1)*novec(2)*novec(3) + kvec(2)*novec(3) + kvec(3))% &
704 interpolator%hierarchy(counter)%parent( &
705 modulo(kd + 1, 2) + stride) = &
706 interpolator%index(lvec(1), lvec(2), lvec(3)) + &
707 kvec(1)*novec(2)*novec(3) + kvec(2)*novec(3) + kvec(3)
709 interpolator%hierarchy( &
710 interpolator%hierarchy(counter)%parent( &
711 modulo(kd + 1, 2) + stride))%children(modulo(kd, 2) + stride) = counter
714 if (interpolator%order == 1)
then
715 interpolator%hierarchy(counter)% &
716 function_type(cdim) = -1
717 elseif (ld == 1 .OR. interpolator%order == 2)
then
718 interpolator%hierarchy(counter)% &
719 function_type(cdim) = 1 + modulo(kd, 2)
720 elseif (ld == 2 .OR. interpolator%order == 3)
then
721 interpolator%hierarchy(counter)% &
722 function_type(cdim) = 3 + modulo(kd, 4)
724 interpolator%hierarchy(counter)% &
725 function_type(cdim) = 7 + modulo(kd, 8)
736 subroutine fg_to_sg(interpolator, fg_values, sg_values)
737 sll_real64,
dimension(:, :, :),
intent(in) :: fg_values
738 sll_real64,
dimension(:),
intent(out) :: sg_values
741 sll_int32,
dimension(3) :: fg_ind
743 do j = 1, interpolator%size_basis
745 sg_values(j) = fg_values(fg_ind(1), fg_ind(2), fg_ind(3));
752 sll_int32,
intent(in) :: sg_index
757 do j = 1, interpolator%dim
758 fg_index(j) = 2**(interpolator%levels(j) - &
759 interpolator%hierarchy(sg_index)%level(j))* &
760 (1 + 2*interpolator%hierarchy(sg_index)%index_on_level(j)) + 1;
773 sll_real64,
dimension(:),
intent(in) :: data_in
774 sll_comp64,
dimension(:),
intent(out) :: data_out
775 sll_int32 :: i1, i2, i3, j
778 do i2 = 0, interpolator%levels(2)
779 do i3 = 0, min(interpolator%max_level - i2, interpolator%levels(3))
780 do j = interpolator%index(0, i2, i3), &
781 interpolator%index(0, i2, i3) + max(2**(i2 - 1), 1)*max(2**(i3 - 1), 1) - 1
782 call interpolator%ToHierarchical1D(1, &
783 min(interpolator%levels(1), interpolator%max_level - i2 - i3), &
784 j, data_in, data_out)
789 do i1 = 0, interpolator%levels(1)
790 do i3 = 0, min(interpolator%max_level - i1, interpolator%levels(3))
791 do j = interpolator%index(i1, 0, i3), &
792 interpolator%index(i1, 0, i3) + max(2**(i1 - 1), 1)*max(2**(i3 - 1), 1) - 1
793 call interpolator%ToHierarchical1D_comp &
795 min(interpolator%levels(2), interpolator%max_level - i1 - i3), j, &
801 do i1 = 0, interpolator%levels(1)
802 do i2 = 0, min(interpolator%max_level - i1, interpolator%levels(2))
803 do j = interpolator%index(i1, i2, 0), &
804 interpolator%index(i1, i2, 0) + max(2**(i1 - 1), 1)*max(2**(i2 - 1), 1) - 1
805 call interpolator%ToHierarchical1D_comp &
807 min(interpolator%levels(3), interpolator%max_level - i1 - i2), j, &
815 subroutine todehi(interpolator, data_array)
817 sll_comp64,
dimension(:),
intent(inout) :: data_array
818 sll_int32 :: i1, i2, i3, j
820 do i1 = 0, interpolator%levels(1)
821 do i2 = 0, min(interpolator%max_level - i1, interpolator%levels(2))
822 do j = interpolator%index(i1, i2, 0), &
823 interpolator%index(i1, i2, 0) + max(2**(i1 - 1), 1)*max(2**(i2 - 1), 1) - 1
824 call interpolator%ToDehi1D &
826 min(interpolator%levels(3), interpolator%max_level - i1 - i2), j, &
832 do i1 = 0, interpolator%levels(1)
833 do i3 = 0, min(interpolator%max_level - i1, interpolator%levels(3))
834 do j = interpolator%index(i1, 0, i3), &
835 interpolator%index(i1, 0, i3) + max(2**(i1 - 1), 1)*max(2**(i3 - 1), 1) - 1
836 call interpolator%ToDehi1D &
838 min(interpolator%levels(2), interpolator%max_level - i1 - i3), j, &
844 do i2 = 0, interpolator%levels(2)
845 do i3 = 0, min(interpolator%max_level - i2, interpolator%levels(3))
846 do j = interpolator%index(0, i2, i3), &
847 interpolator%index(0, i2, i3) + max(2**(i2 - 1), 1)*max(2**(i3 - 1), 1) - 1
848 call interpolator%ToDehi1D(1, &
849 min(interpolator%levels(1), interpolator%max_level - i2 - i3), &
857 subroutine tohira(interpolator, data_array)
859 sll_comp64,
dimension(:),
intent(inout) :: data_array
860 sll_int32 :: i1, i2, i3, j
862 do i2 = 0, interpolator%levels(2)
863 do i3 = 0, min(interpolator%max_level - i2, interpolator%levels(3))
864 do j = interpolator%index(0, i2, i3), &
865 interpolator%index(0, i2, i3) + max(2**(i2 - 1), 1)*max(2**(i3 - 1), 1) - 1
866 call interpolator%ToHira1D(1, &
867 min(interpolator%levels(1), interpolator%max_level - i2 - i3), &
873 do i1 = 0, interpolator%levels(1)
874 do i3 = 0, min(interpolator%max_level - i1, interpolator%levels(3))
875 do j = interpolator%index(i1, 0, i3), &
876 interpolator%index(i1, 0, i3) + max(2**(i1 - 1), 1)*max(2**(i3 - 1), 1) - 1
877 call interpolator%ToHira1D &
879 min(interpolator%levels(2), interpolator%max_level - i1 - i3), j, &
885 do i1 = 0, interpolator%levels(1)
886 do i2 = 0, min(interpolator%max_level - i1, interpolator%levels(2))
887 do j = interpolator%index(i1, i2, 0), &
888 interpolator%index(i1, i2, 0) + max(2**(i1 - 1), 1)*max(2**(i2 - 1), 1) - 1
889 call interpolator%ToHira1D &
891 min(interpolator%levels(3), interpolator%max_level - i1 - i2), j, &
899 subroutine tonodal(interpolator, data_in, data_out)
901 sll_comp64,
dimension(:),
intent(inout) :: data_in
902 sll_real64,
dimension(:),
intent(out) :: data_out
903 sll_int32 :: i1, i2, i3, j
905 do i1 = 0, interpolator%levels(1)
906 do i2 = 0, min(interpolator%max_level - i1, interpolator%levels(2))
907 do j = interpolator%index(i1, i2, 0), &
908 interpolator%index(i1, i2, 0) + max(2**(i1 - 1), 1)*max(2**(i2 - 1), 1) - 1
909 call interpolator%ToNodal1D_comp &
911 min(interpolator%levels(3), interpolator%max_level - i1 - i2), j, &
917 do i1 = 0, interpolator%levels(1)
918 do i3 = 0, min(interpolator%max_level - i1, interpolator%levels(3))
919 do j = interpolator%index(i1, 0, i3), &
920 interpolator%index(i1, 0, i3) + max(2**(i1 - 1), 1)*max(2**(i3 - 1), 1) - 1
921 call interpolator%ToNodal1D_comp &
923 min(interpolator%levels(2), interpolator%max_level - i1 - i3), j, &
929 do i2 = 0, interpolator%levels(2)
930 do i3 = 0, min(interpolator%max_level - i2, interpolator%levels(3))
931 do j = interpolator%index(0, i2, i3), &
932 interpolator%index(0, i2, i3) + max(2**(i2 - 1), 1)*max(2**(i3 - 1), 1) - 1
933 call interpolator%ToNodal1D(1, &
934 min(interpolator%levels(1), interpolator%max_level - i2 - i3), &
935 j, data_in, data_out)
942 subroutine displace(interpolator, dim, displacement, data)
944 sll_comp64,
dimension(:),
intent(inout) ::
data
945 sll_real64,
intent(in) :: displacement
946 sll_int32,
intent(in) :: dim
947 sll_int32 :: i1, i2, i3, j
950 do i2 = 0, interpolator%levels(2)
951 do i3 = 0, min(interpolator%max_level - i2, interpolator%levels(3))
952 do j = interpolator%index(0, i2, i3), &
953 interpolator%index(0, i2, i3) + max(2**(i2 - 1), 1)*max(2**(i3 - 1), 1) - 1
954 call interpolator%Displace1D(1, interpolator%levels(1) - i2 - i3, &
955 j, displacement, data)
959 else if (dim == 2)
then
960 do i1 = 0, interpolator%levels(1)
961 do i3 = 0, min(interpolator%max_level - i1, interpolator%levels(3))
962 do j = interpolator%index(i1, 0, i3), &
963 interpolator%index(i1, 0, i3) + max(2**(i1 - 1), 1)*max(2**(i3 - 1), 1) - 1
964 call interpolator%Displace1D(2, interpolator%levels(2) - i1 - i3, &
965 j, displacement, data)
969 else if (dim == 3)
then
970 do i1 = 0, interpolator%levels(1)
971 do i2 = 0, min(interpolator%max_level - i1, interpolator%levels(2))
972 do j = interpolator%index(i1, i2, 0), &
973 interpolator%index(i1, i2, 0) + max(2**(i1 - 1), 1)*max(2**(i2 - 1), 1) - 1
974 call interpolator%Displace1D(3, interpolator%levels(3) - i1 - i2, &
975 j, displacement, data)
984 subroutine spfft(interpolator, data_in, data_out)
986 sll_real64,
dimension(:),
intent(in) :: data_in
987 sll_comp64,
dimension(:),
intent(out) :: data_out
990 call todehi(interpolator, data_out)
995 subroutine ispfft(interpolator, data_in, data_out)
997 sll_comp64,
dimension(:),
intent(inout) :: data_in
998 sll_real64,
dimension(:),
intent(out) :: data_out
1000 call tohira(interpolator, data_in)
1001 call tonodal(interpolator, data_in, data_out)
Fortran module where set some physical and mathematical constants.
real(kind=f64), parameter, public sll_p_pi
Implementation of a 3D sparse grid with interpolation routines.
subroutine set_hierarchy_info_boundary(interpolator, counter, cdim, lvecin, kvecin, novecin)
Helfer function for initialization. Setting all the information needed for node counter of the sparse...
real(kind=f64) function interpolate_from_hierarchical_surplus_boundary(interpolator, data, eta)
implements interpolation from hierarchical surplus (interpolate_from_interpolant_value) non-periodic
subroutine interpolate_array_disp_sgfft(interpolator, dim, displacment_in, data_in, data_out)
Compute value at displaced grid points using trigonometric interpolation (based on SG FFT)
real(kind=f64) function interpolate_from_hierarchical_surplus(interpolator, data, eta)
Implements interpolate_from_interpolant_value for periodic sparse grid.
subroutine ispfft(interpolator, data_in, data_out)
Sparse grid inverse FFT.
subroutine initialize_sg3d(interpolator, levels, order, interpolation, interpolation_type, eta_min, eta_max, boundary, modified)
Initialization function. Set up the hierarchy of the sparse grid.
real(kind=f64) function interpolate_from_interpolant_value(interpolator, data, eta)
Compute the value of the sparse grid interpolant at position eta (using standard sparse grid interpol...
subroutine interpolate_const_disp(interpolator, dorder, displacement, data_in, data_out, hiera)
Interpolation function for interpolation at (constantly) displaced grid points; displacement only in ...
subroutine tohierarchical(interpolator, data_in, data_out)
subroutine tohira(interpolator, data_array)
integer(kind=i32) function, dimension(3) fg_index(interpolator, sg_index)
Compute the index of a sparse grid node on level "level" with index "index_on_level" on full grid wit...
subroutine fg_to_sg(interpolator, fg_values, sg_values)
Functions to evaluate fg on sg and sg on fg.
subroutine spfft(interpolator, data_in, data_out)
Sparse grid FFT.
subroutine displace(interpolator, dim, displacement, data)
subroutine set_hierarchy_info(interpolator, counter, cdim, lvecin, kvecin, novecin)
Helfer function for initialization. Setting all the information needed for node counter of the sparse...
subroutine todehi(interpolator, data_array)
subroutine tonodal(interpolator, data_in, data_out)
Dimension-independent functions for sparse grid with polynomial basis functions.
Sparse grid object for 3d with interpolation routines.
Class defining the sparse grid data structure.