9 #include "sll_memory.h"
10 #include "sll_working_precision.h"
28 sll_int32,
dimension(:, :),
pointer :: index
62 sll_int32,
intent(in) :: dim
63 sll_real64,
intent(in) :: displacment_in
64 sll_comp64,
dimension(:),
intent(inout) :: data_in
65 sll_real64,
dimension(:),
intent(out) :: data_out
67 sll_real64:: displacement
69 displacement = displacment_in*2.0_f64*
sll_p_pi/interpolator%length(dim)
70 call interpolator%displace(dim, displacement, data_in);
71 call interpolator%ISPFFT(data_in, data_out);
93 sll_real64,
dimension(:),
intent(in) :: eta_min
94 sll_real64,
dimension(:),
intent(in) :: eta_max
95 sll_int32,
dimension(:),
intent(in) :: levels
96 sll_int32,
intent(in) :: order
97 sll_int32,
intent(in) :: interpolation
98 sll_int32,
intent(in) :: interpolation_type
99 sll_int32,
intent(in) :: modified
103 sll_int32,
intent(in) :: boundary
105 sll_int32 :: i, j, k1, k2, l1, l2, l, counter
107 sll_int32,
dimension(:),
allocatable :: novec, lvec, kvec
108 sll_int32 :: no1, no2
110 interpolator%dim = 2;
111 interpolator%modified = modified;
112 interpolator%boundary = boundary;
113 interpolator%max_level = levels(1);
114 do l = 2, interpolator%dim
115 interpolator%max_level = max(levels(l), interpolator%max_level);
117 if (interpolator%modified == 1)
then
118 interpolator%max_level = interpolator%max_level + 1;
121 interpolator%size_basis = 0;
122 if (interpolator%boundary == 0)
then
123 do l = 0, interpolator%max_level
124 do l1 = max(0, l - levels(2)), min(l, levels(1))
126 interpolator%size_basis = &
127 interpolator%size_basis + &
128 max(2**(l1 - 1), 1)*max(2**(l2 - 1), 1);
132 do l = 0, interpolator%max_level
133 do l1 = max(0, l - levels(2)), min(l, levels(1))
145 interpolator%size_basis = &
146 interpolator%size_basis + no1*no2;
151 call interpolator%initialize_sg(levels, order, interpolation, &
152 interpolation_type, eta_min, eta_max);
153 sll_allocate(lvec(interpolator%dim), ierr);
154 sll_allocate(kvec(interpolator%dim), ierr);
155 sll_allocate(novec(interpolator%dim), ierr);
156 sll_allocate(interpolator%index(0:interpolator%levels(1), 0:interpolator%levels(2)), ierr);
159 do l = 0, interpolator%max_level
160 interpolator%level_mapping(l) = counter;
161 do l1 = max(0, l - levels(2)), min(l, levels(1))
162 novec(1) = max(2**(l1 - 1), 1)
163 if (interpolator%boundary == 1 .AND. l1 == 0)
then
168 novec(2) = max(2**(l2 - 1), 1)
169 if (interpolator%boundary == 1 .AND. l2 == 0)
then
173 interpolator%index(l1, l2) = counter
174 do k1 = 0, novec(1) - 1
176 do k2 = 0, novec(2) - 1
178 interpolator%hierarchy(counter)%index_on_level(1) = k1;
179 interpolator%hierarchy(counter)%index_on_level(2) = k2;
180 do j = 1, interpolator%dim
181 if (interpolator%boundary == 0)
then
182 call interpolator%set_hierarchy_info &
183 (counter, j, lvec, kvec, novec);
184 elseif (interpolator%boundary == 1)
then
185 call interpolator%set_hierarchy_info_boundary &
186 (counter, j, lvec, kvec, novec);
189 counter = counter + 1;
194 interpolator%level_mapping(interpolator%max_level + 1) = counter;
196 do i = 1, interpolator%size_basis
197 do j = 1, interpolator%dim
198 interpolator%hierarchy(i)%coordinate(j) = &
199 interpolator%hierarchy(i)%coordinate(j)* &
200 interpolator%length(j) + &
201 interpolator%eta_min(j)
212 sll_int32,
intent(in) :: cdim
213 sll_int32,
intent(in) :: counter
214 sll_int32,
dimension(:),
intent(in) :: lvecin
215 sll_int32,
dimension(:),
intent(in) :: kvecin
216 sll_int32,
dimension(:),
intent(in) :: novecin
217 sll_int32,
dimension(:),
allocatable :: lvec, kvec, novec
218 sll_int32 :: jj, stride
220 sll_allocate(lvec(interpolator%dim), jj);
221 sll_allocate(kvec(interpolator%dim), jj);
222 sll_allocate(novec(interpolator%dim), jj);
223 do jj = 1, interpolator%dim
224 lvec(jj) = lvecin(jj);
225 kvec(jj) = kvecin(jj);
226 novec(jj) = novecin(jj);
230 interpolator%hierarchy(counter)%level(cdim) = ld;
233 interpolator%hierarchy(counter)%coordinate(cdim) = 0.0_f64
234 interpolator%hierarchy(counter)%parent(stride) = &
236 interpolator%hierarchy(counter)%parent(stride + 1) = &
238 interpolator%hierarchy(counter)%function_type(cdim) = 0
240 interpolator%hierarchy(counter)%coordinate(cdim) = &
241 1.0_f64/(2.0_f64**ld) + real(kd, f64)*1.0_f64/(2.0_f64**(ld - 1))
243 lvec(cdim) = lvec(cdim) - 1;
244 novec(cdim) = max(novec(cdim)/2, 1);
246 kvec(cdim) = modulo((kd + 1)/2, max(2**(ld - 2), 1));
247 interpolator%hierarchy(counter)%parent( &
248 modulo(kd, 2) + stride) = &
249 interpolator%hierarchy( &
250 interpolator%index(lvec(1), lvec(2)) + &
251 kvec(1)*novec(2) + kvec(2))%parent(stride)
254 interpolator%hierarchy(counter)%parent( &
255 modulo(kd + 1, 2) + stride) = &
256 interpolator%index(lvec(1), lvec(2)) + &
257 kvec(1)*novec(2) + kvec(2)
259 interpolator%hierarchy( &
260 interpolator%hierarchy(counter)%parent( &
261 modulo(kd + 1, 2) + stride))%children(modulo(kd, 2) + stride) = counter
263 interpolator%hierarchy( &
264 interpolator%hierarchy(counter)%parent( &
265 modulo(kd + 1, 2) + stride))%children(modulo(kd + 1, 2) + stride) = counter
267 if (interpolator%order == 1)
then
268 interpolator%hierarchy(counter)% &
269 function_type(cdim) = -1
270 elseif (ld == 1 .OR. interpolator%order == 2)
then
271 interpolator%hierarchy(counter)% &
272 function_type(cdim) = 1 + modulo(kd, 2)
273 elseif (ld == 2 .OR. interpolator%order == 3)
then
274 interpolator%hierarchy(counter)% &
275 function_type(cdim) = 3 + modulo(kd, 4)
277 interpolator%hierarchy(counter)% &
278 function_type(cdim) = 7 + modulo(kd, 8)
289 sll_int32,
intent(in) :: cdim
290 sll_int32,
intent(in) :: counter
291 sll_int32,
dimension(:),
intent(in) :: lvecin, kvecin, novecin
292 sll_int32,
dimension(:),
allocatable :: lvec, kvec, novec
293 sll_int32 :: jj, stride
295 sll_allocate(lvec(interpolator%dim), jj);
296 sll_allocate(kvec(interpolator%dim), jj);
297 sll_allocate(novec(interpolator%dim), jj);
298 do jj = 1, interpolator%dim
299 lvec(jj) = lvecin(jj);
300 kvec(jj) = kvecin(jj);
301 novec(jj) = novecin(jj);
305 interpolator%hierarchy(counter)%level(cdim) = ld;
308 interpolator%hierarchy(counter)%coordinate(cdim) = real(kd, f64);
309 interpolator%hierarchy(counter)%parent(stride) = &
311 interpolator%hierarchy(counter)%parent(stride + 1) = &
313 interpolator%hierarchy(counter)%function_type(cdim) = -1
315 interpolator%hierarchy(counter)%coordinate(cdim) = &
316 1.0_f64/(2.0_f64**ld) + real(kd, f64)*1.0_f64/(2.0_f64**(ld - 1))
318 lvec(cdim) = lvec(cdim) - 1;
319 novec(cdim) = novec(cdim)/2;
323 interpolator%hierarchy(counter)%parent(stride) = &
324 interpolator%index(lvec(1), lvec(2)) + &
325 kvec(1)*novec(2) + kvec(2)
327 interpolator%hierarchy(counter)%parent( &
328 stride + 1) = interpolator%hierarchy(counter)%parent( &
331 interpolator%hierarchy(counter)%parent( &
332 stride + 1) = interpolator%hierarchy(counter)%parent( &
335 interpolator%hierarchy(interpolator%hierarchy(counter)%parent(stride))%children &
341 kvec(cdim) = modulo((kd + 1)/2, max(2**(ld - 2), 1));
342 kvec(cdim) = (kd + 1)/2;
343 if (kvec(cdim) < novec(cdim))
then
344 interpolator%hierarchy(counter)%parent( &
345 modulo(kd, 2) + stride) = &
346 interpolator%hierarchy( &
347 interpolator%index(lvec(1), lvec(2)) + &
348 kvec(1)*novec(2) + kvec(2))%parent(stride)
350 kvec(cdim) = kvec(cdim) - 1;
351 interpolator%hierarchy(counter)%parent( &
352 modulo(kd, 2) + stride) = &
353 interpolator%hierarchy( &
354 interpolator%index(lvec(1), lvec(2)) + &
355 kvec(1)*novec(2) + kvec(2))%parent(stride + 1)
359 interpolator%hierarchy(counter)%parent( &
360 modulo(kd + 1, 2) + stride) = &
361 interpolator%index(lvec(1), lvec(2)) + &
362 kvec(1)*novec(2) + kvec(2)
364 interpolator%hierarchy( &
365 interpolator%hierarchy(counter)%parent( &
366 modulo(kd + 1, 2) + stride))%children(modulo(kd, 2) + stride) = counter
369 if (interpolator%order == 1)
then
370 interpolator%hierarchy(counter)% &
371 function_type(cdim) = -1
372 elseif (ld == 1 .OR. interpolator%order == 2)
then
373 interpolator%hierarchy(counter)% &
374 function_type(cdim) = 1 + modulo(kd, 2)
375 elseif (ld == 2 .OR. interpolator%order == 3)
then
376 interpolator%hierarchy(counter)% &
377 function_type(cdim) = 3 + modulo(kd, 4)
379 interpolator%hierarchy(counter)% &
380 function_type(cdim) = 7 + modulo(kd, 8)
395 sll_real64,
dimension(:),
intent(in) ::
data
397 sll_real64,
dimension(:),
intent(in) :: eta
399 if (interpolator%boundary == 0)
then
412 sll_int32 :: j, l1, l2, level
414 sll_real64,
dimension(:),
intent(in) ::
data
415 sll_real64,
dimension(:),
intent(in) :: eta
416 sll_real64,
dimension(2) :: eta_norm
417 sll_real64,
dimension(2) :: phi
418 sll_int32,
dimension(2) :: no, l
419 sll_int32,
dimension(:, :),
allocatable :: ind
423 sll_allocate(ind(0:interpolator%max_level, 1:interpolator%dim), j)
426 ind(0:1, 1:interpolator%dim) = 0
428 do j = 1, interpolator%dim
429 eta_norm(j) = (eta(j) - interpolator%eta_min(j))/interpolator%length(j)
430 eta_norm(j) = modulo(eta_norm(j), 1.0_f64)
433 do level = 2, interpolator%levels(j)
434 ind(level, j) = ind(level - 1, j)*2
435 if (eta_norm(j) > scale*real(ind(level, j) + 1, f64))
then
436 ind(level, j) = ind(level, j) + 1
438 scale = scale*0.5_f64
442 do level = 0, interpolator%max_level
443 do l1 = max(0, level - interpolator%levels(2)), min(level, interpolator%levels(1))
445 no(1) = max(2**(l1 - 1), 1)
448 no(2) = max(2**(l2 - 1), 1)
449 index = interpolator%index(l1, l2) + ind(l1, 1)*no(2) &
451 do j = 1, interpolator%dim
452 call interpolator%basis_function(real(2**(max(l(j), 1)), f64)*eta_norm(j) &
453 - real(2*ind(l(j), j), f64) - 1.0_f64, phi(j), &
454 interpolator%hierarchy(index)%function_type(j))
456 val = val +
data(index) &
468 sll_int32 :: j, l1, l2, level
470 sll_real64,
dimension(:),
intent(in) ::
data
471 sll_real64,
dimension(:),
intent(in) :: eta
472 sll_real64,
dimension(2) :: eta_norm
473 sll_real64,
dimension(2) :: phi
474 sll_int32,
dimension(2) :: no, l
475 sll_int32,
dimension(:, :),
allocatable :: ind
476 sll_real64 :: scale, phi1a, phi2a
479 sll_allocate(ind(0:interpolator%max_level, 1:interpolator%dim), j)
482 ind(0:1, 1:interpolator%dim) = 0
484 do j = 1, interpolator%dim
485 eta_norm(j) = (eta(j) - interpolator%eta_min(j))/interpolator%length(j)
488 do level = 2, interpolator%levels(j)
489 ind(level, j) = ind(level - 1, j)*2
490 if (eta_norm(j) > scale*real(ind(level, j) + 1, f64))
then
491 ind(level, j) = ind(level, j) + 1
493 scale = scale*0.5_f64
497 do level = 0, interpolator%max_level
498 do l1 = max(0, level - interpolator%levels(2)), min(level, interpolator%levels(1))
510 no(2) = max(2**(l2 - 1), 1);
514 index = interpolator%index(l1, l2) + ind(l1, 1)*no(2) &
517 do j = 1, interpolator%dim
518 call interpolator%basis_function(real(2**l(j), f64)*eta_norm(j) &
519 - real(2*ind(l(j), j), f64) - 1.0_f64, phi(j), &
520 interpolator%hierarchy(index)%function_type(j))
522 val = val +
data(index) &
525 index = interpolator%index(l1, l2) + ind(l1, 1)*no(2)
526 call interpolator%basis_function(real(2**l(1), f64)*eta_norm(1) &
527 - real(2*ind(l(1), 1), f64) - 1.0_f64, phi(1), &
528 interpolator%hierarchy(index)%function_type(1))
529 call interpolator%basis_function(eta_norm(2), phi(2), -1)
530 val = val +
data(index)*phi(1)*phi(2);
531 call interpolator%basis_function(eta_norm(2) - 1.0_f64, phi(2), -1)
532 val = val +
data(index + 1)*phi(1)*phi(2);
536 index = interpolator%index(l1, l2) + ind(l2, 2)
537 call interpolator%basis_function(real(2**l(2), f64)*eta_norm(2) &
538 - real(2*ind(l(2), 2), f64) - 1.0_f64, phi(2), &
539 interpolator%hierarchy(index)%function_type(2))
540 call interpolator%basis_function(eta_norm(1), phi(1), -1)
541 val = val +
data(index)*phi(1)*phi(2);
542 call interpolator%basis_function(eta_norm(1) - 1.0_f64, phi(1), -1)
543 val = val +
data(index + no(2))*phi(1)*phi(2);
545 index = interpolator%index(l1, l2);
546 call interpolator%basis_function(eta_norm(1), phi(1), -1)
547 call interpolator%basis_function(eta_norm(1) - 1.0_f64, phi1a, -1)
548 call interpolator%basis_function(eta_norm(2), phi(2), -1)
549 call interpolator%basis_function(eta_norm(2) - 1.0_f64, phi2a, -1)
550 val = val +
data(index)*phi(1)*phi(2) + &
551 data(index + 1)*phi(1)*phi2a +
data(index + 2)*phi1a*phi(2) + &
552 data(index + 3)*phi1a*phi2a
567 sll_int32,
dimension(:),
intent(in) :: dorder
568 sll_real64,
dimension(:),
intent(inout) :: data_in
569 sll_real64,
dimension(:),
intent(out) :: data_out
570 sll_real64,
intent(in) ::displacement
571 logical,
intent(in) :: hiera
573 sll_int32 :: j, counter, i2, k2
574 sll_int32,
dimension(2) :: ind_order, no
575 sll_int32,
dimension(:, :),
allocatable :: ind
577 sll_allocate(ind(interpolator%max_level + 1, 2), i2);
578 ind_order(dorder(1)) = 0
580 do i2 = 0, interpolator%levels(dorder(2))
581 ind_order(dorder(2)) = i2
582 no(dorder(2)) = max(2**(i2 - 1), 1);
583 ind(1, dorder(1)) = 0;
584 do k2 = 0, no(dorder(2)) - 1
585 ind(ind_order(dorder(2)) + 1, dorder(2)) = k2;
586 counter = interpolator%index( &
587 ind_order(1), ind_order(2)) + &
588 ind(ind_order(1) + 1, 1)*no(2) &
589 + ind(ind_order(2) + 1, 2)
592 call interpolator%interpolate_disp_1d_periodic &
593 (displacement, dorder(1), &
594 min(interpolator%levels(dorder(1)), interpolator%max_level &
595 - ind_order(dorder(2))), counter, data_in, data_out, hiera);
599 if (hiera .EQV. .false.)
then
601 do j = interpolator%order, 2, -1
602 call interpolator%dehierarchical_part_order &
604 interpolator%dim, 2, dorder, j)
607 call interpolator%dehierarchical_part(data_out, &
608 interpolator%dim, 2, dorder)
620 sll_real64,
dimension(:),
intent(inout) :: data_in
621 sll_real64,
dimension(:),
intent(out) :: data_out
622 sll_int32,
dimension(:),
intent(in) :: dorder
627 sll_real64,
dimension(:),
intent(in) ::displacement
629 sll_int32 :: i1, i2, k2, counter, j, index_parent, index_parent_old
630 sll_int32,
dimension(4) :: no, ind_order
631 sll_int32,
dimension(:, :),
allocatable :: ind
632 sll_real64 :: coordinate_self, coordinate_ancestor, factor, disp
634 sll_allocate(ind(interpolator%max_level + 1, 2), i1);
636 do j = interpolator%order, 2, -1
637 call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part_order &
638 (data_in, 1, 1, dorder, j)
641 call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part(data_in, 1, 1, dorder)
644 ind_order(dorder(1)) = 0
646 do i2 = 0, interpolator%levels(dorder(2))
647 ind_order(dorder(2)) = i2
648 no(dorder(2)) = max(2**(i2 - 1), 1);
649 ind(1, dorder(1)) = 0;
650 do k2 = 0, no(dorder(2)) - 1
651 ind(ind_order(dorder(2)) + 1, dorder(2)) = k2;
652 disp = displacement(2**(i2 - 1) + k2 + 1)
654 counter = interpolator%index( &
655 ind_order(1), ind_order(2)) + &
656 ind(ind_order(1) + 1, 1)*no(2) &
657 + ind(ind_order(2) + 1, 2)
660 call interpolator%sll_t_sparse_grid_interpolator%interpolate_disp_1d_periodic_self &
662 min(interpolator%levels(dorder(1)), interpolator%max_level &
663 - ind_order(dorder(2))), counter, data_in, data_out)
665 index_parent = max( &
666 interpolator%hierarchy(counter)%parent(2*dorder(2) - 1), &
667 interpolator%hierarchy(counter)%parent(2*dorder(2)))
668 coordinate_self = interpolator%hierarchy(counter)%coordinate(dorder(2))
669 index_parent_old = counter;
670 do while (index_parent < index_parent_old)
671 coordinate_ancestor = interpolator%hierarchy(index_parent)% &
672 coordinate(dorder(2))
673 call interpolator%basis_function((coordinate_self - coordinate_ancestor)/ &
674 interpolator%length(dorder(2))* &
675 real(2**(interpolator%hierarchy(index_parent)%level(dorder(2))), f64), &
677 interpolator%hierarchy(index_parent)%function_type(dorder(2)))
678 call interpolator%sll_t_sparse_grid_interpolator%interpolate_disp_1d_periodic_for_neighbor &
680 dorder(1), min(interpolator%levels(dorder(1)), &
681 interpolator%max_level - &
682 interpolator%hierarchy(index_parent)%level(dorder(2))), &
683 min(interpolator%levels(dorder(1)), interpolator%max_level &
684 - ind_order(dorder(2))), index_parent, counter, &
686 index_parent_old = index_parent;
688 max(interpolator%hierarchy(index_parent)%parent(2*dorder(2) - 1), &
689 interpolator%hierarchy(index_parent)%parent(2*dorder(2)))
700 sll_real64,
dimension(:),
intent(inout) :: data_in
701 sll_real64,
dimension(:),
intent(out) :: data_out
702 sll_int32,
dimension(:),
intent(in) :: dorder
703 sll_real64,
intent(in) ::displacement
704 sll_int32 :: i1, i2, k2, counter, j, index_parent, index_parent_old
705 sll_int32,
dimension(4) :: no, ind_order
706 sll_int32,
dimension(:, :),
allocatable :: ind
707 sll_real64 :: coordinate_self, coordinate_ancestor, factor, disp
709 sll_allocate(ind(interpolator%max_level + 1, 2), i1);
711 do j = interpolator%order, 2, -1
712 call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part_order &
713 (data_in, 1, 1, dorder, j)
716 call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part(data_in, 1, 1, dorder)
719 ind_order(dorder(1)) = 0
721 do i2 = 0, interpolator%levels(dorder(2))
722 ind_order(dorder(2)) = i2
723 no(dorder(2)) = max(2**(i2 - 1), 1);
724 ind(1, dorder(1)) = 0;
725 do k2 = 0, no(dorder(2)) - 1
726 ind(ind_order(dorder(2)) + 1, dorder(2)) = k2;
727 counter = interpolator%index( &
728 ind_order(1), ind_order(2)) + &
729 ind(ind_order(1) + 1, 1)*no(2) &
730 + ind(ind_order(2) + 1, 2)
731 disp = displacement*interpolator%hierarchy(counter)%coordinate(dorder(2));
733 call interpolator%sll_t_sparse_grid_interpolator%interpolate_disp_1d_periodic_self &
735 min(interpolator%levels(dorder(1)), interpolator%max_level - &
736 ind_order(dorder(2))), counter, data_in, data_out)
738 index_parent = max( &
739 interpolator%hierarchy(counter)%parent(2*dorder(2) - 1), &
740 interpolator%hierarchy(counter)%parent(2*dorder(2)))
741 coordinate_self = interpolator%hierarchy(counter)%coordinate(dorder(2))
742 index_parent_old = counter;
743 do while (index_parent < index_parent_old)
744 coordinate_ancestor = interpolator%hierarchy(index_parent)% &
745 coordinate(dorder(2))
746 call interpolator%basis_function((coordinate_self - coordinate_ancestor)/ &
747 interpolator%length(dorder(2))* &
748 real(2**(interpolator%hierarchy(index_parent)%level(dorder(2))), f64), &
750 interpolator%hierarchy(index_parent)%function_type(dorder(2)))
751 call interpolator%sll_t_sparse_grid_interpolator%interpolate_disp_1d_periodic_for_neighbor &
753 dorder(1), min(interpolator%levels(dorder(1)), &
754 interpolator%max_level - &
755 interpolator%hierarchy(index_parent)%level(dorder(2))), &
756 min(interpolator%levels(dorder(2)), interpolator%max_level &
757 - ind_order(dorder(2))), index_parent, counter, &
759 index_parent_old = index_parent;
761 max(interpolator%hierarchy(index_parent)%parent(2*dorder(2) - 1), &
762 interpolator%hierarchy(index_parent)%parent(2*dorder(2)))
774 subroutine fg_to_sg(interpolator, fg_values, sg_values)
775 sll_real64,
dimension(:, :),
intent(in) :: fg_values
776 sll_real64,
dimension(:),
intent(out) :: sg_values
780 sll_int32,
dimension(2) :: fg_ind
782 do j = 1, interpolator%size_basis
784 sg_values(j) = fg_values(fg_ind(1), fg_ind(2));
807 sll_int32,
intent(in) :: sg_index
812 do j = 1, interpolator%dim
813 if (interpolator%hierarchy(sg_index)%level(j) == 0)
then
816 fg_index(j) = 2**(interpolator%levels(j) - &
817 interpolator%hierarchy(sg_index)%level(j))* &
818 (1 + 2*interpolator%hierarchy(sg_index)%index_on_level(j)) + 1;
833 sll_real64,
dimension(:),
intent(in) :: data_in
834 sll_comp64,
dimension(:),
intent(out) :: data_out
838 do i = 0, interpolator%levels(2)
839 do j = interpolator%index(0, i), &
840 interpolator%index(0, i) + max(2**(i - 1), 1) - 1
841 call interpolator%ToHierarchical1D(1, &
842 min(interpolator%levels(1), interpolator%max_level - i), &
843 j, data_in, data_out)
846 do i = 0, interpolator%levels(1)
847 do j = interpolator%index(i, 0), &
848 interpolator%index(i, 0) + max(2**(i - 1), 1) - 1
849 call interpolator%ToHierarchical1D_comp(2, &
850 min(interpolator%levels(2), interpolator%max_level - i), j, data_out)
857 subroutine todehi(interpolator, data_array)
859 sll_comp64,
dimension(:),
intent(inout) :: data_array
864 do i = 0, interpolator%levels(1)
865 do j = interpolator%index(i, 0), &
866 interpolator%index(i, 0) + max(2**(i - 1), 1) - 1
867 call interpolator%ToDehi1D(2, &
868 min(interpolator%levels(2), interpolator%max_level - i), j, &
872 do i = 0, interpolator%levels(2)
873 do j = interpolator%index(0, i), &
874 interpolator%index(0, i) + max(2**(i - 1), 1) - 1
875 call interpolator%ToDehi1D(1, &
876 min(interpolator%levels(1), interpolator%max_level - i), j, &
883 subroutine tohira(interpolator, data_array)
885 sll_comp64,
dimension(:),
intent(inout) :: data_array
888 do i = 0, interpolator%levels(2)
889 do j = interpolator%index(0, i), &
890 interpolator%index(0, i) + max(2**(i - 1), 1) - 1
891 call interpolator%ToHira1D(1, &
892 min(interpolator%levels(1), interpolator%max_level - i), j, &
896 do i = 0, interpolator%levels(1)
897 do j = interpolator%index(i, 0), &
898 interpolator%index(i, 0) + max(2**(i - 1), 1) - 1
899 call interpolator%ToHira1D(2, &
900 min(interpolator%levels(2), interpolator%max_level - i), j, &
907 subroutine tonodal(interpolator, data_in, data_out)
909 sll_comp64,
dimension(:),
intent(inout) :: data_in
910 sll_real64,
dimension(:),
intent(out) :: data_out
913 do i = 0, interpolator%levels(1)
914 do j = interpolator%index(i, 0), &
915 interpolator%index(i, 0) + max(2**(i - 1), 1) - 1
916 call interpolator%ToNodal1D_comp(2, &
917 min(interpolator%levels(2), interpolator%max_level - i), &
924 do i = 0, interpolator%levels(2)
925 do j = interpolator%index(0, i), &
926 interpolator%index(0, i) + max(2**(i - 1), 1) - 1
927 call interpolator%ToNodal1D(1, &
928 min(interpolator%levels(1), interpolator%max_level - i), &
929 j, data_in, data_out)
936 subroutine displace(interpolator, dim, displacement, data)
938 sll_comp64,
dimension(:),
intent(inout) ::
data
939 sll_real64,
intent(in) :: displacement
940 sll_int32,
intent(in) :: dim
946 do i = 0, interpolator%levels(2)
947 do j = interpolator%index(0, i), &
948 interpolator%index(0, i) + max(2**(i - 1), 1) - 1
949 call interpolator%Displace1D(1, interpolator%levels(1) - i, &
950 j, displacement, data)
954 do i = 0, interpolator%levels(1)
955 do j = interpolator%index(i, 0), &
956 interpolator%index(i, 0) + max(2**(i - 1), 1) - 1
957 call interpolator%Displace1D(2, interpolator%levels(2) - i, &
958 j, displacement, data)
996 subroutine spfft(interpolator, data_in, data_out)
998 sll_real64,
dimension(:),
intent(in) :: data_in
999 sll_comp64,
dimension(:),
intent(out) :: data_out
1001 call interpolator%ToHierarchical(data_in, data_out)
1002 call interpolator%ToDehi(data_out)
1004 end subroutine spfft
1007 subroutine ispfft(interpolator, data_in, data_out)
1009 sll_comp64,
dimension(:),
intent(inout) :: data_in
1010 sll_real64,
dimension(:),
intent(out) :: data_out
1012 call interpolator%toHira(data_in)
1013 call interpolator%ToNodal(data_in, data_out)
1025 sll_real64,
dimension(:),
intent(inout) ::
data
1026 sll_int32 :: lev, l, j, ind
1028 lev = interpolator%levels(1);
1029 do l = 0, interpolator%max_level - lev
1030 ind = interpolator%index(0, l)
1031 call interpolator%extract_periodic(1, lev, ind,
data, interpolator%stripe)
1032 call interpolator%hierarchical_stripe(interpolator%stripe, lev);
1034 interpolator%stripe(j) = interpolator%stripe(j)/2;
1036 call interpolator%dehierarchical_stripe(interpolator%stripe, lev);
1037 call interpolator%insert_periodic(1, lev, ind, interpolator%stripe, data);
1040 lev = interpolator%levels(2);
1041 do l = 0, interpolator%max_level - lev
1042 ind = interpolator%index(l, 0)
1043 call interpolator%extract_periodic(2, lev, ind,
data, interpolator%stripe)
1044 call interpolator%hierarchical_stripe(interpolator%stripe, lev);
1046 interpolator%stripe(j) = interpolator%stripe(j)*0.5_f64;
1048 call interpolator%dehierarchical_stripe(interpolator%stripe, lev);
1049 call interpolator%insert_periodic(2, lev, ind, interpolator%stripe, data);
1056 sll_real64,
dimension(:),
intent(inout) ::
data
1059 call interpolator%compute_hierarchical_surplus(data);
1060 do j = interpolator%level_mapping(interpolator%max_level) + 1, interpolator%size_basis, 2
1061 data(j) =
data(j)*0.5_f64;
1063 call interpolator%compute_dehierarchical(data);
1082 sll_real64,
dimension(:),
intent(inout) ::
data, hs
1083 sll_real64,
dimension(2),
intent(in) :: width
1084 sll_real64,
dimension(2) :: dx, dxn
1087 do j = 1, interpolator%size_basis
1088 dx = interpolator%hierarchy(j)%coordinate
1089 dxn(1) = dx(1); dxn(2) = dx(2) + width(2)
1090 data(j) =
data(j)*0.5_f64 + 0.125_f64*interpolator%interpolate_from_interpolant_value(hs, dxn);
1091 dxn(2) = dx(2) - width(2)
1092 data(j) =
data(j) + 0.125_f64*interpolator%interpolate_from_interpolant_value(hs, dxn);
1093 dxn(1) = dx(1) + width(1); dxn(2) = dx(2)
1094 data(j) =
data(j) + 0.125_f64*interpolator%interpolate_from_interpolant_value(hs, dxn);
1095 dxn(1) = dx(1) - width(1)
1096 data(j) =
data(j) + 0.125_f64*interpolator%interpolate_from_interpolant_value(hs, dxn);
Fortran module where set some physical and mathematical constants.
real(kind=f64), parameter, public sll_p_pi
Implementation of a 2D sparse grid with interpolation routines.
subroutine displace(interpolator, dim, displacement, data)
Compute the Fourier coefficients of at displaced grid points from Fourier coefficients.
subroutine interpolate_disp_nconst_in_1d(interpolator, displacement, dorder, data_in, data_out)
Functionality: Interpolates the function values for a displacement on in dimension (periodic b....
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 spfft(interpolator, data_in, data_out)
Fourier transform on sparse grid.
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)
subroutine fg_to_sg(interpolator, fg_values, sg_values)
Set sparse grid values from fg vector.
subroutine ispfft(interpolator, data_in, data_out)
Inverse Fourier transform.
subroutine todehi(interpolator, data_array)
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 filter_highest(interpolator, data)
subroutine filter(interpolator, data)
subroutine tohierarchical(interpolator, data_in, data_out)
Compute Fourier coefficient on sparse grid.
subroutine tohira(interpolator, data_array)
integer(kind=i32) function, dimension(2) 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...
real(kind=f64) function interpolate_value_sg(interpolator, data, eta)
Value at eta interpolated from the hierarchical surplus data using standard sparse grid interpolation...
subroutine initialize_sg2d(interpolator, levels, order, interpolation, interpolation_type, eta_min, eta_max, boundary, modified)
Initialize a 2d sparse grid object.
subroutine set_hierarchy_info_boundary(interpolator, counter, cdim, lvecin, kvecin, novecin)
Same as set_hierarchy_info but for points on the boundary along dimension cdim.
real(kind=f64) function interpolate_from_hierarchical_surplus(interpolator, data, eta)
Implementation of interpolate_value_sg for periodic sparse grid.
subroutine tonodal(interpolator, data_in, data_out)
subroutine linear_filter(interpolator, data, hs, width)
subroutine interpolate_disp_linnconst_in_1d(interpolator, displacement, dorder, data_in, data_out)
As interpolate_disp_nconst_in_1d but displacement dependent on displacement*coordinate(dorder(2))
real(kind=f64) function interpolate_from_hierarchical_surplus_boundary(interpolator, data, eta)
Implementation of interpolate_value_sg for sparse grid with boundary.
Dimension-independent functions for sparse grid with polynomial basis functions.
Sparse grid object for 2d with interpolation routines.
Class defining the sparse grid data structure.