Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_sparse_grid_3d.F90
Go to the documentation of this file.
1 
6 
8 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
9 #include "sll_memory.h"
10 #include "sll_working_precision.h"
11 
12  use sll_m_constants, only: &
13  sll_p_pi
14 
17 
18  implicit none
19 
20  public :: &
22 
23  private
24 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
25 
28  sll_int32, dimension(:, :, :), pointer :: index
29 
30  contains
31  procedure :: init => initialize_sg3d! Initialization routine
32  procedure :: interpolate_from_interpolant_value ! Compute the value of the sparse grid interpolant at position eta
34  procedure :: fg_to_sg
35  procedure :: spfft
36  procedure :: ispfft
38 
39  end type sll_t_sparse_grid_3d
40 
41 contains
42 
43 !------------------------------------------------------------------------------!
44 !!!! Interpolation routines !!!!
45 
47  function interpolate_from_interpolant_value(interpolator, data, eta) result(val)
48  class(sll_t_sparse_grid_3d), intent(inout) :: interpolator
49  sll_real64 :: val
50  sll_real64, dimension(:), intent(in) :: data
51  sll_real64, dimension(:), intent(in) :: eta
52 
53  if (interpolator%boundary == 0) then
54  val = interpolate_from_hierarchical_surplus(interpolator, data, eta)
55  else
56  val = interpolate_from_hierarchical_surplus_boundary(interpolator, data, eta)
57  end if
59 
65  subroutine interpolate_const_disp(interpolator, dorder, displacement, data_in, data_out, hiera)
66  class(sll_t_sparse_grid_3d), intent(inout) :: interpolator
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
72 
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
76 
77  sll_allocate(ind(interpolator%max_level + 1, 3), i1);
78  ind_order(dorder(1)) = 0
79  no(dorder(1)) = 1
80  do i3 = 0, interpolator%levels(dorder(3))
81  ind_order(dorder(3)) = i3
82  l(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)
97  ! Evaluate along dorder(1)-stripe
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)
103  end do
104  end do
105  end do
106  end do
107 
108  if (hiera .EQV. .false.) then
109  ! Dehierarchization along dimension dorder(1) only
110  do j = interpolator%order, 2, -1
111  call interpolator%dehierarchical_part_order &
112  (data_out, &
113  interpolator%dim, 2, dorder, j)
114  end do
115 
116  call interpolator%dehierarchical_part(data_out, &
117  interpolator%dim, 2, dorder)
118  end if
119 
120  end subroutine interpolate_const_disp
121 
122 ! helper functions
124  function interpolate_from_hierarchical_surplus(interpolator, data, eta) result(val)
125  class(sll_t_sparse_grid_3d), intent(inout) :: interpolator
126  sll_int32 :: j, l1, l2, l3, level
127  sll_real64 :: val
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
133  sll_real64 :: scale
134  sll_int32 :: index
135 
136  sll_allocate(ind(0:interpolator%max_level, 1:3), j)
137 
138  val = 0.0_f64
139  ind(0:1, 1:3) = 0
140 
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)
144 
145  scale = 0.5_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
150  end if
151  scale = scale*0.5_f64
152  end do
153  end do
154 
155  do level = 0, interpolator%max_level
156  do l1 = 0, min(level, interpolator%levels(1))
157  l(1) = l1
158  no(1) = max(2**(l1 - 1), 1)
159  do l2 = max(0, level - l1 - interpolator%levels(3)), min(level - l1, interpolator%levels(2))
160  l(2) = l2
161  no(2) = max(2**(l2 - 1), 1)
162  l(3) = level - l1 - l2
163  l3 = l(3);
164  no(3) = max(2**(l(3) - 1), 1)
165 
166  index = interpolator%index(l1, l2, l3) + ind(l1, 1)*no(2)*no(3) &
167  + ind(l2, 2)*no(3) + ind(l3, 3)
168  do j = 1, 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))
172  end do
173  val = val + data(index)*phi(1)*phi(2)*phi(3)
174  end do
175  end do
176  end do
177 
179 
181  function interpolate_from_hierarchical_surplus_boundary(interpolator, data, eta) result(val)
182  class(sll_t_sparse_grid_3d), intent(inout) :: interpolator
183  sll_int32 :: j, l1, l2, l3, level
184  sll_real64 :: val
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
192  sll_int32 :: index
193 
194  sll_allocate(ind(0:interpolator%max_level, 1:interpolator%dim), j)
195 
196  val = 0.0_f64
197  ind(0:1, 1:interpolator%dim) = 0
198 
199  do j = 1, interpolator%dim
200  eta_norm(j) = (eta(j) - interpolator%eta_min(j))/interpolator%length(j)
201 
202  scale = 0.5_f64
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
207  end if
208  scale = scale*0.5_f64
209  end do
210  end do
211 
212  do level = 0, interpolator%max_level
213  do l1 = 0, min(level, interpolator%levels(1))
214  l(1) = l1
215  if (l1 == 0) then
216  no(1) = 2;
217  else
218  no(1) = 2**(l1 - 1);
219  end if
220  do l2 = max(0, level - l1 - interpolator%levels(3)), min(level - l1, interpolator%levels(2))
221  l(2) = l2
222  if (l2 == 0) then
223  no(2) = 2;
224  else
225  no(2) = max(2**(l2 - 1), 1);
226  end if
227  l3 = level - l1 - l2;
228  if (l3 == 0) then
229  no(3) = 2;
230  else
231  no(3) = max(2**(l3 - 1), 1);
232  end if
233  if (l1 > 0) then
234  if (l2 > 0) then
235  if (l3 > 0) then
236  index = interpolator%index(l1, l2, l3) + &
237  ind(l1, 1)*no(2)*no(3) &
238  + ind(l2, 2)*no(3) + ind(l3, 3)
239 
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))
244  end do
245  val = val + data(index) &
246  *phi(1)*phi(2)*phi(3)
247  else ! l3=0
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);
261  end if
262  else !l2 = 0
263  if (l3 > 0) then
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);
277  else !l3=l2=0
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
291  end if
292  end if
293  else ! l1=0
294  if (l2 > 0) then
295  if (l3 > 0) then
296  index = interpolator%index(l1, l2, l3) + &
297  ind(l1, 1)*no(2)*no(3) &
298  + ind(l2, 2)*no(3) + ind(l3, 3)
299 
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))
304  end do
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)
309  else ! l1=l3=0
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
324  end if
325  else !l1=l2 = 0!AB hier noch nicht ganz ueberarbeitet
326  if (l3 > 0) then
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);
341  else !l3=l2=0
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
358  end if
359  end if
360  end if
361  end do
362  end do
363  end do
365 
366 !!!! End interpolation routines !!!!
367 !------------------------------------------------------------------------------!
368 
369 !------------------------------------------------------------------------------!
370 !!!! SGFFT routines !!!!!
371 
373  subroutine interpolate_array_disp_sgfft(interpolator, dim, displacment_in, data_in, data_out)
374  class(sll_t_sparse_grid_3d), intent(inout) :: interpolator
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
380 
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);
384  end subroutine interpolate_array_disp_sgfft
385 
386 !!!! End SGFFT routines !!!!!
387 !------------------------------------------------------------------------------!
388 
389 !------------------------------------------------------------------------------!
390 !!!! Initialization routines !!!!
391 
393  subroutine initialize_sg3d( &
394  interpolator, &
395  levels, &
396  order, &
397  interpolation, &
398  interpolation_type, &
399  eta_min, &
400  eta_max, &
401  boundary, &
402  modified)
403 
404  class(sll_t_sparse_grid_3d), intent(inout) :: interpolator
405  sll_real64, dimension(:), intent(in) :: eta_min
406  sll_real64, dimension(:), intent(in) :: eta_max
407  sll_int32, dimension(:), intent(in) :: levels
408 
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
413 
417  sll_int32, intent(in) :: boundary
418 
419 
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
423 
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);
433  end do
434  if (interpolator%modified == 1) then
435  interpolator%max_level = interpolator%max_level + 2;
436  end if
437 
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))
443  l3 = l - l1 - l2
444  interpolator%size_basis = &
445  interpolator%size_basis + &
446  max(2**(l1 - 1), 1)*max(2**(l2 - 1), 1)*max(2**(l3 - 1), 1)
447  end do
448  end do
449  end do
450  else
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))
454  l3 = l - l1 - l2
455  if (l1 == 0) then
456  no1 = 2;
457  else
458  no1 = 2**(l1 - 1);
459  end if
460  if (l2 == 0) then
461  no2 = 2;
462  else
463  no2 = 2**(l2 - 1);
464  end if
465  if (l3 == 0) then
466  no3 = 2;
467  else
468  no3 = 2**(l3 - 1);
469  end if
470  interpolator%size_basis = &
471  interpolator%size_basis + no1*no2*no3;
472  end do
473  end do
474  end do
475  end if
476 
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)
480 
481  ! Set the hierarchy of the grid
482  counter = 1
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
488  novec(1) = 2;
489  end if
490  lvec(1) = l1;
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
494  novec(2) = 2;
495  end if
496  lvec(2) = l2;
497  lvec(3) = l - l1 - l2;
498  l3 = lvec(3);
499  novec(3) = max(2**(l3 - 1), 1)
500  if (interpolator%boundary == 1 .AND. l3 == 0) then
501  novec(3) = 2;
502  end if
503  lvec(3) = l3
504  interpolator%index(l1, l2, l3) = counter
505  do k1 = 0, novec(1) - 1
506  kvec(1) = k1
507  do k2 = 0, novec(2) - 1
508  kvec(2) = k2
509  do k3 = 0, novec(3) - 1
510  kvec(3) = k3
511  do j = 1, interpolator%dim
512  if (interpolator%boundary == 0) then
513  call set_hierarchy_info(interpolator, counter, j, &
514  lvec, kvec, novec);
515  else
517  (interpolator, counter, j, lvec, kvec, novec);
518  end if
519  end do
520  counter = counter + 1;
521  end do
522  end do
523  end do
524  end do
525  end do
526  end do
527  interpolator%level_mapping(interpolator%max_level + 1) = counter;
528  ! Now rescale all the coordinates to the actual mesh size
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)
535  end do
536  end do
537 
538  end subroutine initialize_sg3d
539 
542  subroutine set_hierarchy_info(interpolator, counter, cdim, lvecin, kvecin, novecin)
543  class(sll_t_sparse_grid_3d), intent(inout) :: interpolator
544  sll_int32 :: ld
545  sll_int32 :: kd
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
551 
552  sll_int32, dimension(3) :: lvec, kvec, novec
553  sll_int32 :: jj, stride
554 
555  do jj = 1, interpolator%dim
556  lvec(jj) = lvecin(jj);
557  kvec(jj) = kvecin(jj);
558  novec(jj) = novecin(jj);
559  end do
560  ld = lvec(cdim);
561  kd = kvec(cdim);
562  interpolator%hierarchy(counter)%level(cdim) = ld;
563  interpolator%hierarchy(counter)%index_on_level(cdim) = kd;
564  stride = cdim*2 - 1
565  if (ld == 0) then
566  interpolator%hierarchy(counter)%coordinate(cdim) = 0.0_f64
567  interpolator%hierarchy(counter)%parent(stride) = &
568  counter
569  interpolator%hierarchy(counter)%parent(stride + 1) = &
570  counter
571  interpolator%hierarchy(counter)%function_type(cdim) = 0
572  else
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))
575 
576  lvec(cdim) = lvec(cdim) - 1;
577  novec(cdim) = max(novec(cdim)/2, 1);
578  ! This one is actually only a neighboring point not the direct parent.
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) + &
585  kvec(2)*novec(3) + &
586  kvec(3))%parent(stride)
587  ! This is the actual parent.
588  kvec(cdim) = kd/2;
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) + &
593  kvec(3)
594  ! Now tell my parent that I am his child.
595  interpolator%hierarchy( &
596  interpolator%hierarchy(counter)%parent( &
597  modulo(kd + 1, 2) + stride))%children(modulo(kd, 2) + stride) = counter
598  if (ld == 1) then
599  interpolator%hierarchy( &
600  interpolator%hierarchy(counter)%parent( &
601  modulo(kd + 1, 2) + stride))%children(modulo(kd + 1, 2) + stride) = counter
602  end if
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 !(order==3) then!
610  interpolator%hierarchy(counter)% &
611  function_type(cdim) = 3 + modulo(kd, 4)
612  else
613  interpolator%hierarchy(counter)% &
614  function_type(cdim) = 7 + modulo(kd, 8)
615  end if
616  end if
617 
618  end subroutine set_hierarchy_info
619 
621  subroutine set_hierarchy_info_boundary(interpolator, counter, cdim, lvecin, kvecin, novecin)
622  class(sll_t_sparse_grid_3d), intent(inout) :: interpolator
623  sll_int32 :: ld
624  sll_int32 :: kd
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
630 
631  sll_int32, dimension(:), allocatable :: lvec, kvec, novec
632  sll_int32 :: jj, stride
633 
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);
641  end do
642  ld = lvec(cdim);
643  kd = kvec(cdim);
644  interpolator%hierarchy(counter)%level(cdim) = ld;
645  stride = cdim*2 - 1
646  if (ld == 0) then
647  interpolator%hierarchy(counter)%coordinate(cdim) = real(kd, f64);
648  interpolator%hierarchy(counter)%parent(stride) = &
649  counter
650  interpolator%hierarchy(counter)%parent(stride + 1) = &
651  counter
652  interpolator%hierarchy(counter)%function_type(cdim) = -1
653  else
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))
656 
657  lvec(cdim) = lvec(cdim) - 1;
658  novec(cdim) = novec(cdim)/2;
659  if (ld == 1) then
660  kvec(cdim) = kd/2;
661  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)
665  if (cdim == 1) then
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( &
672  stride) + novec(3)
673  else
674  interpolator%hierarchy(counter)%parent( &
675  stride + 1) = interpolator%hierarchy(counter)%parent( &
676  stride) + 1
677  end if
678  interpolator%hierarchy(interpolator%hierarchy(counter)%parent(stride))%children &
679  (stride) = counter
680  !interpolator%hierarchy(interpolator%hierarchy(counter)%parent(stride+1))%&
681  ! children(stride) = counter
682  else
683  ! This one is actually only a neighboring point not the direct parent.
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))% &
692  parent(stride)
693  else
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))% &
700  parent(stride + 1)
701  end if
702  ! This is the actual parent.
703  kvec(cdim) = kd/2;
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)
708  ! Now tell my parent that I am his child.
709  interpolator%hierarchy( &
710  interpolator%hierarchy(counter)%parent( &
711  modulo(kd + 1, 2) + stride))%children(modulo(kd, 2) + stride) = counter
712  end if
713 
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 !(order==3) then!
721  interpolator%hierarchy(counter)% &
722  function_type(cdim) = 3 + modulo(kd, 4)
723  else
724  interpolator%hierarchy(counter)% &
725  function_type(cdim) = 7 + modulo(kd, 8)
726  end if
727  end if
728 
729  end subroutine set_hierarchy_info_boundary
730 
731 !!!! End initialization routines !!!!
732 !------------------------------------------------------------------------------!
733 
734 !------------------------------------------------------------------------------!
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
739  class(sll_t_sparse_grid_3d), intent(in) :: interpolator
740  sll_int32 :: j
741  sll_int32, dimension(3) :: fg_ind
742 
743  do j = 1, interpolator%size_basis
744  fg_ind = fg_index(interpolator, j);
745  sg_values(j) = fg_values(fg_ind(1), fg_ind(2), fg_ind(3));
746  end do
747 
748  end subroutine fg_to_sg
749 
751  function fg_index(interpolator, sg_index)
752  sll_int32, intent(in) :: sg_index
753  sll_int32, dimension(3) :: fg_index
754  class(sll_t_sparse_grid_3d), intent(in) :: interpolator
755  sll_int32 :: j
756 
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;
761  end do
762 
763  end function fg_index
764 
765 ! End functions fg_to_sg and sg_to_fg
766 !------------------------------------------------------------------------------!
767 
768 !------------------------------------------------------------------------------!
769 !!!! SGFFT helper functions (Some clean-up needed) !!!!
770 
771  subroutine tohierarchical(interpolator, data_in, data_out)
772  class(sll_t_sparse_grid_3d), intent(inout) :: interpolator
773  sll_real64, dimension(:), intent(in) :: data_in
774  sll_comp64, dimension(:), intent(out) :: data_out
775  sll_int32 :: i1, i2, i3, j
776 
777 ! Maybe possible with index in index to make dimension independent
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)
785  end do
786  end do
787  end do
788 
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 &
794  (2, &
795  min(interpolator%levels(2), interpolator%max_level - i1 - i3), j, &
796  data_out)
797  end do
798  end do
799  end do
800 
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 &
806  (3, &
807  min(interpolator%levels(3), interpolator%max_level - i1 - i2), j, &
808  data_out)
809  end do
810  end do
811  end do
812 
813  end subroutine tohierarchical
814 
815  subroutine todehi(interpolator, data_array)
816  class(sll_t_sparse_grid_3d), intent(inout) :: interpolator
817  sll_comp64, dimension(:), intent(inout) :: data_array
818  sll_int32 :: i1, i2, i3, j
819 
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 &
825  (3, &
826  min(interpolator%levels(3), interpolator%max_level - i1 - i2), j, &
827  data_array)
828  end do
829  end do
830  end do
831 
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 &
837  (2, &
838  min(interpolator%levels(2), interpolator%max_level - i1 - i3), j, &
839  data_array)
840  end do
841  end do
842  end do
843 
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), &
850  j, data_array)
851  end do
852  end do
853  end do
854 
855  end subroutine todehi
856 
857  subroutine tohira(interpolator, data_array)
858  class(sll_t_sparse_grid_3d), intent(inout) :: interpolator
859  sll_comp64, dimension(:), intent(inout) :: data_array
860  sll_int32 :: i1, i2, i3, j
861 
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), &
868  j, data_array)
869  end do
870  end do
871  end do
872 
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 &
878  (2, &
879  min(interpolator%levels(2), interpolator%max_level - i1 - i3), j, &
880  data_array)
881  end do
882  end do
883  end do
884 
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 &
890  (3, &
891  min(interpolator%levels(3), interpolator%max_level - i1 - i2), j, &
892  data_array)
893  end do
894  end do
895  end do
896 
897  end subroutine tohira
898 
899  subroutine tonodal(interpolator, data_in, data_out)
900  class(sll_t_sparse_grid_3d), intent(inout) :: interpolator
901  sll_comp64, dimension(:), intent(inout) :: data_in
902  sll_real64, dimension(:), intent(out) :: data_out
903  sll_int32 :: i1, i2, i3, j
904 
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 &
910  (3, &
911  min(interpolator%levels(3), interpolator%max_level - i1 - i2), j, &
912  data_in)
913  end do
914  end do
915  end do
916 
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 &
922  (2, &
923  min(interpolator%levels(2), interpolator%max_level - i1 - i3), j, &
924  data_in)
925  end do
926  end do
927  end do
928 
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)
936  end do
937  end do
938  end do
939 
940  end subroutine tonodal
941 
942  subroutine displace(interpolator, dim, displacement, data)
943  class(sll_t_sparse_grid_3d), intent(inout) :: interpolator
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
948 
949  if (dim == 1) then
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)
956  end do
957  end do
958  end do
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)
966  end do
967  end do
968  end do
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)
976  end do
977  end do
978  end do
979  end if
980 
981  end subroutine displace
982 
984  subroutine spfft(interpolator, data_in, data_out)
985  class(sll_t_sparse_grid_3d), intent(inout) :: interpolator
986  sll_real64, dimension(:), intent(in) :: data_in
987  sll_comp64, dimension(:), intent(out) :: data_out
988 
989  call tohierarchical(interpolator, data_in, data_out)
990  call todehi(interpolator, data_out)
991 
992  end subroutine spfft
993 
995  subroutine ispfft(interpolator, data_in, data_out)
996  class(sll_t_sparse_grid_3d), intent(inout) :: interpolator
997  sll_comp64, dimension(:), intent(inout) :: data_in
998  sll_real64, dimension(:), intent(out) :: data_out
999 
1000  call tohira(interpolator, data_in)
1001  call tonodal(interpolator, data_in, data_out)
1002 
1003  end subroutine ispfft
1004 
1005 !!!! End SGFFT helper functions !!!!
1006 !------------------------------------------------------------------------------!
1007 
1008 end module sll_m_sparse_grid_3d
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.
    Report Typos and Errors