Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_sparse_grid_interpolator.F90
Go to the documentation of this file.
1 
6 module sll_m_sparse_grid_interpolator
7 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8 #include "sll_memory.h"
9 #include "sll_working_precision.h"
10 #include "sll_fftw.h"
11 
12  use iso_c_binding, only: &
13  c_double_complex, &
14  c_size_t
15 
16  use sll_m_fft, only: &
17  sll_t_fft, &
18  sll_s_fft_init_c2c_1d, &
19  sll_s_fft_exec_c2c_1d, &
20  sll_s_fft_free, &
21  sll_f_fft_allocate_aligned_complex, &
22  sll_p_fft_forward, &
23  sll_p_fft_backward, &
24  sll_p_fft_measure
25 
26  use sll_m_interpolators_1d_base, only: &
27  sll_c_interpolator_1d
28 
29  use sll_m_periodic_interpolator_1d, only: &
30  sll_t_periodic_interpolator_1d
31 
32  implicit none
33 
34  public :: &
35  sll_t_sparse_grid_interpolator
36 
37  private
38 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
39 
41  type :: fft_hierarchical
42  type(sll_t_fft) :: fw
43  type(sll_t_fft) :: bw
44  complex(c_double_complex), dimension(:), pointer :: in
45  complex(c_double_complex), dimension(:), pointer :: out
46  integer(c_size_t) :: sz_in
47  integer(c_size_t) :: sz_out
48  end type fft_hierarchical
49 
51  type sparsegrid_node
52  sll_real64, dimension(:), allocatable :: coordinate
53  sll_int32, dimension(:), allocatable :: parent
54  ! type(sparsegrid_node_ptr),dimension(:), pointer :: parents
55  sll_int32, dimension(:), allocatable :: children
56  sll_int32, dimension(:), allocatable :: function_type
57  sll_int32, dimension(:), allocatable :: level
58  sll_int32, dimension(:), allocatable :: index_on_level
59  !sll_real64, dimension(:), allocatable :: values
60  end type sparsegrid_node
61 
62  type :: interpolator_base_ptr
63  class(sll_c_interpolator_1d), pointer :: ptr
64  end type interpolator_base_ptr
65 
67  type :: sll_t_sparse_grid_interpolator
68  sll_int32 :: dim
69  sll_real64, dimension(:), pointer :: eta_min
70  sll_real64, dimension(:), pointer :: eta_max
71  sll_real64, dimension(:), pointer :: length
72  sll_real64 :: volume
73  sll_int32 :: max_level
74  sll_int32, dimension(:), pointer :: levels
75  sll_int32 :: order
76  sll_int32 :: interpolation
77  sll_int32 :: no_basis_functions
78  sll_int32 :: size_basis
79  sll_int32 :: modified
80  sll_int32 :: boundary
81 
82  class(sparsegrid_node), dimension(:), pointer :: hierarchy
83 
84  class(fft_hierarchical), dimension(:), pointer :: fft_object
85 
86  !sll_int32, dimension(:,:,:,:), pointer :: index
87  sll_real64, dimension(:), pointer :: stripe, stripe_out
88  sll_real64, dimension(:), pointer :: hs_weights
89  sll_int32, dimension(:), pointer :: hs_weights_index
90  type(sll_t_periodic_interpolator_1d), dimension(:, :), pointer :: interp_per
91  !real,dimension(:), target :: interp_per_x
92  !type(per_1d_interpolator),dimension(:,:), pointer :: interp_v
93  ! type(odd_degree_spline_1d_interpolator),dimension(:,:), pointer :: interp_v
94  !type(lagrange_1d_interpolator),dimension(:,:), pointer :: interpl_v
95  type(interpolator_base_ptr), dimension(:, :), pointer :: interp
96  !type(sll_c_interpolator_1d), dimension(:,:), pointer :: interp
97  sll_int32, dimension(:), pointer :: level_mapping
98 
99  contains
100  procedure :: compute_hierarchical_surplus
101  procedure :: compute_linear_hierarchical_surplus
102  procedure :: compute_dehierarchical
103  procedure :: dehierarchical
104  procedure :: dehierarchical_part
105  procedure :: hierarchical_part
106  procedure :: interpolate_disp
107  procedure :: integrate_trapezoidal
108  procedure :: integrate_trapezoidal2
109  procedure :: extract_periodic
110  procedure :: hierarchical_stripe
111  procedure :: dehierarchical_stripe
112  procedure :: insert_periodic
113  procedure :: interpolate_disp_1d_periodic_self
114  procedure, nopass :: basis_function
115  procedure, nopass :: basis_function_derivative
116  procedure :: dehierarchical_part_order
117  procedure :: interpolate_disp_1d_periodic_for_neighbor
118  procedure :: interpolate_disp_1d_periodic
119  procedure :: displace1d
120  procedure :: tonodal1d
121  procedure :: tonodal1d_comp
122  procedure :: todehi1d
123  procedure :: tohira1d
124  procedure :: tohierarchical1d
125  procedure :: tohierarchical1d_comp
126  procedure :: initialize_sg
127  procedure :: hierarchical_part_order
128  procedure :: dehierarchical_order
129  procedure :: free => free_sparse_grid
130 
131  end type sll_t_sparse_grid_interpolator
132 
133 contains
134 
135 !------------------------------------------------------------------------------!
136 !!!! Helper functions for initialization
137 
138  subroutine initialize_sg( &
139  interpolator, &
140  levels, &
141  order, &
142  interpolation, &
143  interpolation_type, &
144  eta_min, &
145  eta_max)
146  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
147 
148  sll_real64, dimension(:), intent(in) :: eta_min
149  sll_real64, dimension(:), intent(in) :: eta_max
150  sll_int32, dimension(:), intent(in) :: levels
151  sll_int32, intent(in) :: order
152  sll_int32, intent(in) :: interpolation
153  sll_int32, intent(in) :: interpolation_type
154 
155  sll_int32 :: i, j
156  sll_int32 :: ierr
157 
158  sll_allocate(interpolator%eta_min(interpolator%dim), ierr);
159  sll_allocate(interpolator%eta_max(interpolator%dim), ierr);
160  sll_allocate(interpolator%length(interpolator%dim), ierr);
161  sll_allocate(interpolator%levels(interpolator%dim), ierr);
162  !SLL_ALLOCATE( interpolator, ierr )
163  interpolator%volume = 1.0_f64;
164  do j = 1, interpolator%dim
165  interpolator%eta_min(j) = eta_min(j);
166  interpolator%eta_max(j) = eta_max(j);
167  interpolator%length(j) = eta_max(j) - eta_min(j);
168  interpolator%volume = interpolator%volume*interpolator%length(j);
169  end do
170  interpolator%levels = levels;
171  sll_allocate(interpolator%level_mapping(0:interpolator%max_level + 1), ierr);
172  interpolator%order = order
173  interpolator%interpolation = interpolation;
174  if (order == 1) then
175  interpolator%no_basis_functions = 2
176  elseif (order == 2) then
177  interpolator%no_basis_functions = 3
178  elseif (order == 3) then
179  interpolator%no_basis_functions = 7
180  elseif (order == 4) then
181  interpolator%no_basis_functions = 15
182  end if
183 
184  print *, "Size of the basis: ", interpolator%size_basis
185  sll_allocate(interpolator%hierarchy(interpolator%size_basis), ierr)
186  do j = 1, interpolator%size_basis
187  sll_allocate(interpolator%hierarchy(j)%coordinate(interpolator%dim), ierr);
188  sll_allocate(interpolator%hierarchy(j)%parent(2*interpolator%dim), ierr);
189  sll_allocate(interpolator%hierarchy(j)%children(2*interpolator%dim), ierr);
190  sll_allocate(interpolator%hierarchy(j)%level(interpolator%dim), ierr);
191  sll_allocate(interpolator%hierarchy(j)%index_on_level(interpolator%dim), ierr);
192  sll_allocate(interpolator%hierarchy(j)%function_type(interpolator%dim), ierr);
193  end do
194 
195  sll_allocate(interpolator%stripe(2**interpolator%max_level + 1), ierr)
196  sll_allocate(interpolator%stripe_out(2**interpolator%max_level + 1), ierr)
197 
198  do j = 1, interpolator%size_basis
199  interpolator%hierarchy(j)%children = -1
200  end do
201 
202  if (interpolator%order > 1) then
203  sll_allocate(interpolator%hs_weights_index(interpolator%order + 1), ierr)
204  interpolator%hs_weights_index(1) = 1
205  interpolator%hs_weights_index(2) = 1
206  interpolator%hs_weights_index(3) = 3
207  if (interpolator%order == 2) then
208  sll_allocate(interpolator%hs_weights(2), ierr)
209  elseif (interpolator%order == 3) then
210  interpolator%hs_weights_index(4) = 7
211  sll_allocate(interpolator%hs_weights(6), ierr)
212  interpolator%hs_weights(3) = -0.125_f64
213  interpolator%hs_weights(4) = -interpolator%hs_weights(3)
214  interpolator%hs_weights(5) = interpolator%hs_weights(4)
215  interpolator%hs_weights(6) = interpolator%hs_weights(3)
216  elseif (interpolator%order == 4) then
217  interpolator%hs_weights_index(4) = 7
218  interpolator%hs_weights_index(5) = 15
219  sll_allocate(interpolator%hs_weights(14), ierr)
220  interpolator%hs_weights(3) = -0.125_f64
221  interpolator%hs_weights(4) = -interpolator%hs_weights(3)
222  interpolator%hs_weights(5) = interpolator%hs_weights(4)
223  interpolator%hs_weights(6) = interpolator%hs_weights(3)
224  interpolator%hs_weights(7) = -1.0_f64/16.0_f64
225  interpolator%hs_weights(8) = 5.0_f64/112.0_f64
226  interpolator%hs_weights(9) = interpolator%hs_weights(7)
227  interpolator%hs_weights(10) = 7.0_f64/80.0_f64
228  interpolator%hs_weights(11) = 7.0_f64/80.0_f64
229  interpolator%hs_weights(12) = interpolator%hs_weights(7)
230  interpolator%hs_weights(13) = interpolator%hs_weights(8)
231  interpolator%hs_weights(14) = interpolator%hs_weights(7)
232  end if
233  interpolator%hs_weights(1) = -0.25_f64
234  interpolator%hs_weights(2) = interpolator%hs_weights(1)
235  end if
236 
237  sll_allocate(interpolator%interp_per(interpolator%dim, interpolator%max_level), ierr);
238  !SLL_ALLOCATE(interpolator%interp_v(interpolator%dim/2,interpolator%max_level),ierr);
239  !SLL_ALLOCATE(interpolator%interpl_v(interpolator%dim/2,interpolator%max_level),ierr);
240  sll_allocate(interpolator%interp(interpolator%dim, interpolator%max_level), ierr);
241  do i = 1, interpolator%max_level
242  do j = 1, interpolator%dim
243  if (interpolation_type == 0) then
244  call interpolator%interp_per(j, i)%init(2**i + 1, &
245  interpolator%eta_min(j), &
246  interpolator%eta_max(j), &
247  1, interpolation);
248  ! call interpolator%interp_v(j,i)%init( 2**i + 1, &
249  ! interpolator%eta_min(j+1), &
250  ! interpolator%eta_max(j+1),&
251  ! interpolation-1);
252 
253  else
254  call interpolator%interp_per(j, i)%init(2**i + 1, &
255  interpolator%eta_min(j), &
256  interpolator%eta_max(j), &
257  2, interpolation);
258  !call interpolator%interpl_v(j,i)%init( 2**i +1,&
259  ! interpolator%eta_min(j+interpolator%dim/2), &
260  ! interpolator%eta_max(j+interpolator%dim/2),&
261  ! HERMITE_LAGRANGE,interpolation/2);
262  end if
263 
264  ! interpolator%interp(j,i)=interpolator%interp_x(j,i);
265  ! interpolator%interp(j+interpolator%dim/2,i) = interpolator%interp_v(j,i);
266  end do
267  end do
268 
269  ! For FFT
270  ALLOCATE (interpolator%fft_object(interpolator%max_level + 1))
271  call fft_initialize(interpolator%fft_object, interpolator%max_level)
272 
273  end subroutine initialize_sg
274 
275 !!!! End helper functions for initialization
276 !------------------------------------------------------------------------------!
277 
278 !------------------------------------------------------------------------------!
279 !!!! Functions to build hierarchical surplus !!!!
280 
281 ! Public function to be called
282  subroutine compute_hierarchical_surplus(interpolator, data_array)
283  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
284  sll_real64, dimension(:), intent(inout) :: data_array
285  sll_int32 :: i
286 
287  call hierarchical(interpolator, data_array);
288  do i = 2, interpolator%order
289  call hierarchical_order(interpolator, data_array, i);
290  end do
291 
292  end subroutine compute_hierarchical_surplus
293 
294  subroutine compute_linear_hierarchical_surplus(interpolator, data_array)
295  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
296  sll_real64, dimension(:), intent(inout) :: data_array
297 
298  call hierarchical(interpolator, data_array);
299  end subroutine compute_linear_hierarchical_surplus
300 
301 ! Public function to be called
302  subroutine compute_dehierarchical(interpolator, data_array)
303  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
304  sll_real64, dimension(:), intent(inout) :: data_array
305  sll_int32 :: i
306 
307  do i = interpolator%order, 2, -1
308  call dehierarchical_order(interpolator, data_array, i);
309  end do
310 
311  call dehierarchical(interpolator, data_array);
312  end subroutine compute_dehierarchical
313 
314 !!!! End functions to build hierarchical surplus !!!!
315 !------------------------------------------------------------------------------!
316 
317 !------------------------------------------------------------------------------!
318 !!!! Helper functions for interpolation routines !!!!
319 
320 !!!! End helper functions for interpolation routines !!!!
321 !------------------------------------------------------------------------------!
322 
323 !------------------------------------------------------------------------------!
324 !!!! Evaluation of basis functions !!!!
325 
326  subroutine basis_function(x, fx, type)
327  sll_real64, intent(in) :: x
328  sll_real64, intent(inout) :: fx
329  sll_int32, intent(in) :: type
330 
331  select case (type)
332  case (-1)
333  !print*, -1
334  if ((x >= -1.0_f64) .AND. (x < 0.0_f64)) then
335  fx = 1.0_f64 + x
336  elseif ((x >= 0.0_f64) .AND. (x < 1.0_f64)) then
337  fx = 1.0_f64 - x
338  else
339  fx = 0.0_f64
340  end if
341  case (0)
342  !print*, 0
343  fx = 1.0_f64; ! CHANGE_CONSTANT
344  !if ((x>=-1.0_f64) .AND. (x<0.0_f64)) then
345  ! fx = -x
346  !elseif ((x>=0.0_f64) .AND. (x<1.0_f64)) then
347  ! fx = x
348  !else
349  ! fx = 0.0_f64
350  !end if
351  case (1, 2)
352  if ((x > -1.0_f64) .AND. (x < 1.0_f64)) then
353  fx = (1.0_f64 - x)*(x + 1.0_f64)
354  else
355  fx = 0.0_f64
356  end if
357  case (3, 5)
358  if ((x > -1.0_f64) .AND. (x < 1.0_f64)) then
359  fx = (1.0_f64 - x)*(x + 1.0_f64)*(1.0_f64 - x/3.0_f64)
360  else
361  fx = 0.0_f64
362  end if
363  case (4, 6)
364  if ((x > -1.0_f64) .AND. (x < 1.0_f64)) then
365  fx = (1.0_f64 - x)*(x + 1.0_f64)*(1.0_f64 + x/3.0_f64)
366  else
367  fx = 0.0_f64
368  end if
369  case (7, 11)
370  if ((x > -1.0_f64) .AND. (x < 1.0_f64)) then
371  fx = (1.0_f64 - x)*(1.0_f64 + x)*(1.0_f64 - x/3.0_f64)*(1.0_f64 - x/7.0_f64)
372  else
373  fx = 0.0_f64
374  end if
375  case (8, 12)
376  if ((x > -1.0_f64) .AND. (x < 1.0_f64)) then
377  fx = (1.0_f64 - x)*(1.0_f64 + x)*(1.0_f64 + x/3.0_f64)*(1.0_f64 - x/5.0_f64)
378  else
379  fx = 0.0_f64
380  end if
381  case (9, 13)
382  if ((x > -1.0_f64) .AND. (x < 1.0_f64)) then
383  fx = (1.0_f64 - x)*(1.0_f64 + x)*(1.0_f64 - x/3.0_f64)*(1.0_f64 + x/5.0_f64)
384  else
385  fx = 0.0_f64
386  end if
387  case (10, 14)
388  if ((x > -1.0_f64) .AND. (x < 1.0_f64)) then
389  fx = (1.0_f64 - x)*(1.0_f64 + x)*(1.0_f64 + x/3.0_f64)*(1.0_f64 + x/7.0_f64)
390  else
391  fx = 0.0_f64
392  end if
393  end select
394  end subroutine basis_function
395 
396  subroutine basis_function_derivative(x, fx, type)
397  sll_real64, intent(in) :: x
398  sll_real64, intent(inout) :: fx
399  sll_int32, intent(in) :: type
400 
401  select case (type)
402  case (0)
403  fx = 0.0_f64! CHANGE_CONSTANT
404  !if ((x>=-1.0_f64) .AND. (x<0.0_f64)) then
405  ! fx = -1.0_f64
406  !elseif ((x>=0.0_f64) .AND. (x<1.0_f64)) then
407  ! fx = 1.0_f64
408  !else
409  ! fx = 0.0_f64
410  !end if
411  case (1)
412  if ((x >= -1.0_f64) .AND. (x < 1.0_f64)) then
413  fx = -2.0_f64*x
414  else
415  fx = 0.0_f64
416  end if
417  case (2)
418  if ((x >= -1.0_f64) .AND. (x < 1.0_f64)) then
419  fx = -1.0_f64/3.0_f64 - 2.0_f64*x + x*x
420  else
421  fx = 0.0_f64
422  end if
423  case (3)
424  if ((x >= -1.0_f64) .AND. (x < 1.0_f64)) then
425  fx = 1.0_f64/3.0_f64 - 2.0_f64*x - x*x
426  else
427  fx = 0.0_f64
428  end if
429  end select
430 
431  end subroutine basis_function_derivative
432 
433 !!!! End evaluation of basis functions !!!!
434 !------------------------------------------------------------------------------!
435 
436 !------------------------------------------------------------------------------!
437 !!!! Helper functions on 1D stripes (extract, insert, displace, interpolate) !!!
438 
439 ! Displace functions for 1D
440  subroutine displace_on_stripe_periodic_for_neighbor(interpolator, displacement, dim, max_level, max_level_neighbor)
441  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
442  sll_real64, intent(in) :: displacement
443  sll_int32, intent(in) :: max_level, max_level_neighbor, dim
444  sll_int32 :: cell, size_fraction, size_neighbor
445 
446  call displace_on_stripe_periodic(interpolator, displacement, dim, max_level)
447 
448  size_neighbor = 2**max_level_neighbor
449  size_fraction = 2**max_level/size_neighbor
450 
451  do cell = 1, size_neighbor
452  interpolator%stripe_out(cell) = interpolator%stripe_out(1 + (cell - 1)*size_fraction)
453  end do
454 
455  end subroutine displace_on_stripe_periodic_for_neighbor
456 
457 ! end functions for neighbors
458 
459  subroutine displace_on_stripe_periodic(interpolator, displacement, dim, max_level)
460  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
461  sll_real64, intent(in) :: displacement
462  sll_int32, intent(in) :: max_level, dim
463  sll_int32 :: size
464 
465  size = 2**max_level + 1;
466  if (max_level == 0) then
467  interpolator%stripe_out(1) = interpolator%stripe(1)
468  else
469  interpolator%stripe(size) = interpolator%stripe(1);
470  call interpolator%interp_per(dim,max_level)%interpolate_array_disp(size, interpolator%stripe(1:size), displacement, interpolator%stripe_out(1:size))
471  end if
472 
473  end subroutine displace_on_stripe_periodic
474 
475 ! Interpolator functions in 1D
476 
477 ! Interpolate along dimension dim only (periodic boundary conditions), interpolation at a displaced function value
478  subroutine interpolate_disp_1d_periodic(interpolator, displacement, dim, max_level, index, data_in, data_out, hiera)
479  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
480  sll_real64, dimension(:), intent(in) :: data_in
481  sll_real64, dimension(:), intent(out) :: data_out
482  sll_real64, intent(in) :: displacement
483  sll_int32, intent(in) :: dim, max_level, index
484  logical, intent(in) :: hiera
485 
486  call extract_periodic(interpolator, dim, max_level, &
487  index, data_in, interpolator%stripe)
488 
489  call dehierarchical_stripe(interpolator, &
490  interpolator%stripe, max_level)
491 
492  call displace_on_stripe_periodic(interpolator, displacement, dim, max_level)
493 
494  if (hiera) then
495  call hierarchical_stripe(interpolator, &
496  interpolator%stripe_out, max_level);
497  end if
498 
499  call insert_periodic(interpolator, dim, max_level, &
500  index, interpolator%stripe_out, data_out)
501 
502  end subroutine interpolate_disp_1d_periodic
503 
504  recursive subroutine interpolate_disp_recursive(interpolator, no_dims, dim, node_index, displacement, data_in, data_out, hiera)
505  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
506  sll_real64, intent(in) :: displacement
507  sll_real64, dimension(:), intent(in) :: data_in
508  sll_real64, dimension(:), intent(out) :: data_out
509  sll_int32, intent(in) :: dim, no_dims
510  sll_int32, intent(in) :: node_index
511  logical, intent(in) :: hiera
512  sll_int32 :: j, child_index, max_level
513 
514  max_level = interpolator%max_level!interpolator%levels(dim)
515  do j = 1, no_dims
516  if (j .NE. dim) then
517  max_level = max_level - interpolator%hierarchy(node_index)%level(j);
518  end if
519  end do
520  max_level = min(max_level, interpolator%levels(dim))
521 
522  call interpolate_disp_1d_periodic(interpolator, displacement, dim, &
523  max_level, node_index, data_in, data_out, hiera)
524 
525  do j = 1, 2*no_dims
526  if ((j + 1)/2 .NE. dim) then
527  child_index = interpolator%hierarchy(node_index)%children(j);
528  if (child_index > 0) then
529  call interpolate_disp_recursive(interpolator, no_dims, dim, &
530  child_index, displacement, data_in, data_out, hiera)
531  end if
532  end if
533  end do
534 
535  end subroutine interpolate_disp_recursive
536 
537 ! Interpolation function for interpolation at (constantly) displaced grid points; displacement only in dimension dim
538  subroutine interpolate_disp(interpolator, dim, displacement, data_in, data_out, hiera)
539  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
540  sll_real64, dimension(:), intent(inout) :: data_in
541  sll_real64, dimension(:), intent(out) :: data_out
542  sll_int32, intent(in) :: dim
543  sll_real64, intent(in) ::displacement
544  sll_int32 :: j, jj
545  sll_int32, dimension(:), allocatable :: dorder
546  logical, intent(in) :: hiera
547 
548  sll_allocate(dorder(interpolator%dim), j);
549  dorder(1) = dim;
550  jj = 2;
551  do j = 1, interpolator%dim
552  if (j .NE. dim) then
553  dorder(jj) = j;
554  jj = jj + 1;
555  end if
556  end do
557 
558  call interpolate_disp_recursive(interpolator, interpolator%dim, dim, 1, &
559  displacement, data_in, data_out, hiera);
560  if (hiera .EQV. .false.) then
561  ! Dehierarchization along dimension dorder(1) only
562  do j = interpolator%order, 2, -1
563  call dehierarchical_part_order &
564  (interpolator, data_out, &
565  interpolator%dim, 2, dorder, j)
566  end do
567 
568  call dehierarchical_part(interpolator, data_out, &
569  interpolator%dim, 2, dorder)
570  end if
571 
572 !!$ SLL_ALLOCATE(dorder(interpolator%dim),j);
573 !!$ dorder = (1:interpolator%dim);
574 !!$ dorder(1) = dim;
575 !!$ jj = 2;
576 !!$ do j=1,interpolator%dim
577 !!$ if(j .NE. dim) then
578 !!$ dorder(jj) = j;
579 !!$ jj = jj+1;
580 !!$ end if
581 !!$ end do
582 !!$
583 !!$ call interpolate_disp_recursive_self(interpolator,interpolator%dim,dim,1,displacement, data_in, data_out)
584 !!$
585 !!$
586 !!$ ! Dehierarchization along dimension dorder(1) only
587 !!$ do j=interpolator%order,2,-1
588 !!$ call dehierarchical_part_order&
589 !!$ (interpolator%sparse_grid_interpolator,data_in,&
590 !!$ interpolator%dim,2,dorder,j)
591 !!$ end do
592 !!$
593 !!$ call dehierarchical_part(interpolator%sparse_grid_interpolator,data_in,&
594 !!$ interpolator%dim,2,dorder)
595 
596  end subroutine interpolate_disp
597 
598 !!!!!
599 ! Interpolates along dimensions "dim" on the stripe defined by "index". Then result is only computed for every "size_fraction" element and inserted into the global vector for the stripe defined by "index_neighbor".
600 !!!!!!!!
601 subroutine interpolate_disp_1d_periodic_for_neighbor(interpolator,displacement,factor,dim,max_level,max_level_neighbor,index,index_neighbor,data_in,data_out)
602  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
603  sll_real64, dimension(:), intent(in) :: data_in
604  sll_real64, dimension(:), intent(out) :: data_out
605  sll_real64, intent(in) :: displacement, factor
606  sll_int32, intent(in) :: dim, max_level, index, index_neighbor, max_level_neighbor
607 
608  call extract_periodic(interpolator, dim, max_level, &
609  index, data_in, interpolator%stripe)
610  !call dehi_periodic(interpolator,&
611  ! interpolator%stripe,max_level)
612 
613  call displace_on_stripe_periodic_for_neighbor(interpolator, displacement, dim, &
614  max_level, max_level_neighbor)
615 
616  !if(index_neighbor == 4) then
617  ! print*, interpolator%stripe_out(1:2**max_level_neighbor)
618  !end if
619 
620  !call hira_periodic (interpolator,&
621  ! interpolator%stripe_out, max_level)
622  call insert_periodic_additive(interpolator, factor, dim, max_level_neighbor, &
623  index_neighbor, interpolator%stripe_out, data_out)
624 
625  end subroutine interpolate_disp_1d_periodic_for_neighbor
626 
627  subroutine interpolate_disp_1d_periodic_self(interpolator, displacement, dim, max_level, index, data_in, data_out)
628  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
629  sll_real64, dimension(:), intent(in) :: data_in
630  sll_real64, dimension(:), intent(out) :: data_out
631  sll_real64, intent(in) :: displacement
632  sll_int32, intent(in) :: dim, max_level, index
633 
634  call extract_periodic(interpolator, dim, max_level, &
635  index, data_in, interpolator%stripe)
636 
637  call displace_on_stripe_periodic(interpolator, displacement, dim, max_level)
638 
639  !call hira_periodic (interpolator,&
640  ! interpolator%stripe_out, max_level)
641  call insert_periodic(interpolator, dim, max_level, &
642  index, interpolator%stripe_out, data_out)
643 
644  end subroutine interpolate_disp_1d_periodic_self
645 
646 ! Extract functions (extract 1D stripe from dD)
647 
648  subroutine extract_periodic(sparsegrid, dim, max_level, index, data_in, data_out)
649  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
650  sll_int32, intent(in) :: dim, max_level, index
651  sll_int32 :: n_points, index_running
652  sll_real64, dimension(:), intent(in) :: data_in
653  sll_real64, dimension(:), intent(out) :: data_out
654 
655  n_points = 2**(max_level)
656  data_out(1) = data_in(index)
657  if (max_level > 0) then
658  index_running = sparsegrid%hierarchy(index)%children(dim*2)
659  data_out(n_points/2 + 1) = data_in(index_running)
660  if (max_level > 1) then
661  call extract_recursive(sparsegrid, n_points/2 + 1, n_points/4, &
662  index_running, dim, data_in, data_out)
663  end if
664  end if
665 
666  end subroutine extract_periodic
667 
668  recursive subroutine extract_recursive(sparsegrid, index_stripe, stride, index_sg, dim, data_in, data_out)
669  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
670  sll_int32, intent(in) :: index_stripe, stride, index_sg, dim
671  sll_real64, dimension(:), intent(in) :: data_in
672  sll_real64, dimension(:), intent(inout) :: data_out
673  !print*, 'expr', index_sg,index_stripe
674 
675  data_out(index_stripe - stride) = &
676  data_in(sparsegrid%hierarchy(index_sg)%children(dim*2 - 1))
677  data_out(index_stripe + stride) = &
678  data_in(sparsegrid%hierarchy(index_sg)%children(dim*2))
679  if (stride > 1) then
680  call extract_recursive(sparsegrid, index_stripe - stride, stride/2, &
681  sparsegrid%hierarchy(index_sg)%children(dim*2 - 1), dim, &
682  data_in, data_out)
683  call extract_recursive(sparsegrid, index_stripe + stride, stride/2, &
684  sparsegrid%hierarchy(index_sg)%children(dim*2), dim, &
685  data_in, data_out)
686  end if
687  end subroutine extract_recursive
688 
689 ! Insert functions (write 1D stripe back to dD)
690 
691  subroutine insert_periodic(sparsegrid, dim, max_level, index, data_in, data_out)
692  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
693  sll_int32, intent(in) :: dim, max_level, index
694  sll_int32 :: n_points, index_running
695  sll_real64, dimension(:), intent(in) :: data_in
696  sll_real64, dimension(:), intent(inout) :: data_out
697 
698  n_points = 2**(max_level)
699  data_out(index) = data_in(1)
700  if (max_level > 0) then
701  index_running = sparsegrid%hierarchy(index)%children(dim*2)
702 
703  data_out(index_running) = data_in(n_points/2 + 1)
704  if (max_level > 1) then
705 
706  call insert_recursive(sparsegrid, n_points/2 + 1, n_points/4, index_running, &
707  dim, data_in, data_out)
708  end if
709  end if
710 
711  end subroutine insert_periodic
712 
713  recursive subroutine insert_recursive(sparsegrid, index_stripe, stride, index_sg, dim, data_in, data_out)
714  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
715  sll_int32, intent(in) :: index_stripe, stride, index_sg, dim
716  sll_real64, dimension(:), intent(in) :: data_in
717  sll_real64, dimension(:), intent(inout) :: data_out
718 
719  data_out(sparsegrid%hierarchy(index_sg)%children(dim*2 - 1)) = &
720  (data_in(index_stripe - stride))
721  data_out(sparsegrid%hierarchy(index_sg)%children(dim*2)) = &
722  (data_in(index_stripe + stride))
723  if (stride > 1) then
724  call insert_recursive(sparsegrid, index_stripe - stride, stride/2, &
725  sparsegrid%hierarchy(index_sg)%children(dim*2 - 1), dim, &
726  data_in, data_out)
727  call insert_recursive(sparsegrid, index_stripe + stride, stride/2, &
728  sparsegrid%hierarchy(index_sg)%children(dim*2), dim, &
729  data_in, data_out)
730  end if
731  end subroutine insert_recursive
732 
733  subroutine insert_periodic_additive(sparsegrid, factor, dim, max_level, index, data_in, data_out)
734  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
735  sll_real64, intent(in) :: factor
736  sll_int32, intent(in) :: dim, max_level, index
737  sll_int32 :: n_points, index_running
738  sll_real64, dimension(:), intent(in) :: data_in
739  sll_real64, dimension(:), intent(inout) :: data_out
740 
741  n_points = 2**(max_level)
742  data_out(index) = data_out(index) + data_in(1)*factor
743  if (max_level > 0) then
744  index_running = sparsegrid%hierarchy(index)%children(dim*2)
745 
746  data_out(index_running) = data_out(index_running) + data_in(n_points/2 + 1)*factor
747  if (max_level > 1) then
748 
749  call insert_additive_recursive(sparsegrid, factor, n_points/2 + 1, n_points/4, index_running, &
750  dim, data_in, data_out)
751  end if
752  end if
753 
754  end subroutine insert_periodic_additive
755 
756  recursive subroutine insert_additive_recursive(sparsegrid, factor, index_stripe, stride, index_sg, dim, data_in, data_out)
757  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
758  sll_real64, intent(in) :: factor
759  sll_int32, intent(in) :: index_stripe, stride, index_sg, dim
760  sll_real64, dimension(:), intent(in) :: data_in
761  sll_real64, dimension(:), intent(inout) :: data_out
762 
763  data_out(sparsegrid%hierarchy(index_sg)%children(dim*2 - 1)) = &
764  data_out(sparsegrid%hierarchy(index_sg)%children(dim*2 - 1)) + (data_in(index_stripe - stride))*factor
765  !TODO: Should this be there or not?
766  !if(sparsegrid%hierarchy(index_sg)%children(dim*2) .NE. sparsegrid%hierarchy(index_sg)%children(dim*2-1)) then
767  data_out(sparsegrid%hierarchy(index_sg)%children(dim*2)) = &
768  data_out(sparsegrid%hierarchy(index_sg)%children(dim*2)) + (data_in(index_stripe + stride))*factor
769  !end if
770  if (stride > 1) then
771  call insert_additive_recursive(sparsegrid, factor, index_stripe - stride, stride/2, &
772  sparsegrid%hierarchy(index_sg)%children(dim*2 - 1), dim, &
773  data_in, data_out)
774  call insert_additive_recursive(sparsegrid, factor, index_stripe + stride, stride/2, &
775  sparsegrid%hierarchy(index_sg)%children(dim*2), dim, &
776  data_in, data_out)
777  end if
778  end subroutine insert_additive_recursive
779 
780 !!!! End helper function on 1D stripes !!!!
781 !------------------------------------------------------------------------------!
782 
783 !------------------------------------------------------------------------------!
784 !!!! General sparse grid helper functions !!!!
785 
786 ! Various dehierarchical routines (for computing hierarchical surplus)
787 
788 ! Hierarchization of the sparse grid linear hierarchical surplus (data on out) (along all dimensions)
789  subroutine hierarchical(interpolator, data)
790  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
791  sll_real64, dimension(:), intent(inout) :: data
792  !sll_real64, dimension(:), intent(out) :: data_out
793  sll_int32 :: counter
794  sll_real64 :: factor
795 
796  do counter = interpolator%size_basis, 1, -1
797  !write(16,*), 'Now:', l
798  !data_out(counter) = 0.0_f64!data_array(counter,1)
799  factor = 1.0_f64
800  call dehierarchical_d_dimension &
801  (interpolator, data(counter), &
802  data, interpolator%hierarchy(counter)%level, factor, counter, 0)
803  end do
804  end subroutine hierarchical
805 
806 ! Hierarchization of the higher order sparse grid hierarchical surplus (data on output) (along all dimensions)
807  subroutine hierarchical_order(interpolator, data, order)
808  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
809  sll_real64, dimension(:), intent(inout) :: data
810  sll_int32, intent(in) :: order
811  sll_int32 :: counter, start_level, order_level_size, j
812  sll_int32, dimension(:), allocatable :: k
813  sll_real64 :: factor
814 
815  start_level = order - 1
816  order_level_size = 2**(order - 1)
817  sll_allocate(k(interpolator%dim), j);
818  do counter = interpolator%size_basis, 1, -1
819  factor = 1.0_f64
820 
821  do j = 1, interpolator%dim
822  k(j) = modulo(interpolator%hierarchy(counter)%function_type(j) - &
823  interpolator%hs_weights_index(order), &
824  order_level_size) + interpolator%hs_weights_index(order);
825  end do
826  call dehierarchical_order_d_dimension &
827  (interpolator, data(counter), &
828  data, start_level, &
829  interpolator%hierarchy(counter)%level, &
830  k, &
831  factor, counter, 0)
832  end do
833 
834  end subroutine hierarchical_order
835 
836 ! Dehierarchization of the sparse grid linear hierarchical surplus (data) (along all dimensions)
837  subroutine dehierarchical(interpolator, data)
838  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
839  sll_real64, dimension(:), intent(inout) :: data
840  !sll_real64, dimension(:), intent(out) :: data_out
841  sll_int32 :: counter
842  sll_real64 :: factor
843 
844  do counter = 1, interpolator%size_basis
845  !write(16,*), 'Now:', l
846  !data_out(counter) = 0.0_f64!data_array(counter,1)
847  factor = -1.0_f64
848  call dehierarchical_d_dimension &
849  (interpolator, data(counter), &
850  data, interpolator%hierarchy(counter)%level, factor, counter, 0)
851  end do
852  end subroutine dehierarchical
853 
854 ! Dehierarchization of the higher order sparse grid hierarchical surplus (data) (along all dimensions)
855  subroutine dehierarchical_order(interpolator, data, order)
856  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
857  sll_real64, dimension(:), intent(inout) :: data
858  sll_int32, intent(in) :: order
859  sll_int32 :: counter, start_level, order_level_size, j
860  sll_int32, dimension(:), allocatable :: k
861  sll_real64 :: factor
862 
863  start_level = order - 1
864  order_level_size = 2**(order - 1)
865  sll_allocate(k(interpolator%dim), j);
866  do counter = 1, interpolator%size_basis
867  factor = -1.0_f64
868 
869  do j = 1, interpolator%dim
870  k(j) = modulo(interpolator%hierarchy(counter)%function_type(j) - &
871  interpolator%hs_weights_index(order), &
872  order_level_size) + interpolator%hs_weights_index(order);
873  end do
874  call dehierarchical_order_d_dimension &
875  (interpolator, data(counter), &
876  data, start_level, &
877  interpolator%hierarchy(counter)%level, &
878  k, &
879  factor, counter, 0)
880  end do
881 
882  end subroutine dehierarchical_order
883 
884 ! Recursive worker function for dehierarchical
885  recursive subroutine dehierarchical_d_dimension(interpolator, surplus, data_array, level, factor, index, d)
886  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
887  sll_real64, intent(inout) :: surplus
888  sll_real64, dimension(:), intent(in) :: data_array
889  sll_int32, dimension(:), intent(in) :: level
890  sll_real64, intent(in) :: factor
891  sll_int32, intent(in) :: index
892  sll_int32, intent(in) :: d
893  sll_int32 :: j
894 
895  if (index .NE. -1) then
896 
897  if (d > 0) then
898  surplus = surplus + data_array(index)*factor
899  end if
900  do j = d + 1, interpolator%dim
901  if (level(j) > 0) then
902  call dehierarchical_d_dimension(interpolator, surplus, data_array, &
903  level, &
904  -0.5_f64*factor, interpolator%hierarchy(index)%parent(2*j - 1), j)
905  call dehierarchical_d_dimension(interpolator, surplus, data_array, &
906  level, &
907  -0.5_f64*factor, interpolator%hierarchy(index)%parent(2*j), j)
908  end if
909  end do
910  end if
911 
912  end subroutine dehierarchical_d_dimension
913 
914 ! Recursive worker function for dehierarchical order
915  recursive subroutine dehierarchical_order_d_dimension(interpolator, surplus, data_array, start_level, level, k, factor, index, d)
916  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
917  sll_real64, intent(inout) :: surplus
918  sll_real64, dimension(:), intent(in) :: data_array
919  sll_int32, intent(in) :: start_level
920  !sll_real64,dimension(:), intent(in) :: weights
921  sll_int32, dimension(:), intent(in) :: level, k
922  sll_real64, intent(in) :: factor
923  sll_int32, intent(in) :: index
924  sll_int32, intent(in) :: d
925  sll_int32 :: j, father
926 
927  if (d > 0) then
928  surplus = surplus + data_array(index)*factor
929  end if
930  do j = d + 1, interpolator%dim
931  if (level(j) > start_level) then
932  father = max(interpolator%hierarchy(index)%parent(2*j - 1), &
933  interpolator%hierarchy(index)%parent(2*j))
934  if (father .NE. -1) then
935  call dehierarchical_order_d_dimension( &
936  interpolator, surplus, data_array, start_level, level, k, &
937  interpolator%hs_weights(k(j))*factor, father, j)
938  end if
939  end if
940  end do
941 
942  end subroutine dehierarchical_order_d_dimension
943 
944 !!! (De)hierarchization functions on a 1D stripe
945 
946  subroutine hierarchical_stripe(sparsegrid, data, max_level)
947  class(sll_t_sparse_grid_interpolator), intent(inout) :: sparsegrid
948  sll_real64, dimension(:), intent(inout) :: data
949  sll_int32, intent(in) :: max_level
950  sll_int32 :: index, stride, index_run, j, upper, od, level, weights_index, weights_number, factor
951 
952  stride = 1;
953  index = 2;
954  upper = 2**max_level;
955  do level = max_level, 1, -1
956  index_run = index
957  do j = 1, 2**(level - 1)
958  data(index_run) = data(index_run) - &
959  data(index_run - stride)*0.5_f64 - &
960  data(modulo(index_run + stride - 1, upper) + 1)*0.5_f64
961  index_run = index_run + 2*stride
962  end do
963  index = index + stride
964  stride = stride*2
965  end do
966 
967  do od = 2, sparsegrid%order
968  weights_index = sparsegrid%hs_weights_index(od);
969  weights_number = sparsegrid%hs_weights_index(od + 1) - &
970  sparsegrid%hs_weights_index(od)
971  stride = 1;
972  index = 2;
973  do level = max_level, od, -1
974  index_run = index
975  factor = 1
976  do j = 1, 2**(level - 1)
977  data(index_run) = data(index_run) + &
978  data(index_run + factor*stride)* &
979  sparsegrid%hs_weights(weights_index + modulo(j - 1, weights_number))
980  index_run = index_run + 2*stride
981  factor = -factor
982  end do
983  index = index + stride
984  stride = stride*2
985  end do
986  end do
987 
988  end subroutine hierarchical_stripe
989 
990  subroutine dehierarchical_stripe(sparsegrid, data, max_level)
991  class(sll_t_sparse_grid_interpolator), intent(inout) :: sparsegrid
992  sll_real64, dimension(:), intent(inout) :: data
993  sll_int32, intent(in) :: max_level
994  sll_int32 :: index_stripe, stride, od, weights_index, weights_number, index, factor, level, j, index_run
995 
996  do od = sparsegrid%order, 2, -1
997  weights_index = sparsegrid%hs_weights_index(od);
998  weights_number = sparsegrid%hs_weights_index(od + 1) - sparsegrid%hs_weights_index(od)
999  stride = 2**(max_level - od);
1000  index = stride + 1;
1001  do level = od, max_level
1002  index_run = index
1003  factor = 1
1004  do j = 1, 2**(level - 1)
1005  data(index_run) = data(index_run) - &
1006  data(index_run + factor*stride)* &
1007  sparsegrid%hs_weights(weights_index + modulo(j - 1, weights_number))
1008  index_run = index_run + 2*stride
1009  factor = -factor
1010  end do
1011  stride = stride/2
1012  index = index - stride
1013  end do
1014  end do
1015 
1016  if (max_level > 0) then
1017  index_stripe = 2**(max_level - 1)
1018 
1019  data(index_stripe + 1) = data(index_stripe + 1) + data(1)
1020  if (max_level > 1) then
1021  index_stripe = index_stripe/2
1022 
1023  call dehierarchical_stripe_recursive(index_stripe + 1, index_stripe, index_stripe*4, data)
1024 
1025  call dehierarchical_stripe_recursive(index_stripe*3 + 1, index_stripe, index_stripe*4, data)
1026  end if
1027  end if
1028  end subroutine dehierarchical_stripe
1029 
1030  recursive subroutine dehierarchical_stripe_order_recursive(index, stride, data_out)
1031  sll_int32, intent(in) :: index, stride
1032  sll_real64, dimension(:), intent(inout) :: data_out
1033 
1034  data_out(index - stride) = data_out(index - stride) + 0.25_f64*data_out(index)
1035  data_out(index + stride) = data_out(index + stride) + 0.25_f64*data_out(index)
1036 
1037  if (stride > 1) then
1038  call dehierarchical_stripe_order_recursive(index - stride, stride/2, data_out)
1039  call dehierarchical_stripe_order_recursive(index + stride, stride/2, data_out)
1040  end if
1041 
1042  end subroutine dehierarchical_stripe_order_recursive
1043 
1044  recursive subroutine dehierarchical_stripe_recursive(index, stride, upper, data_out)
1045  sll_int32, intent(in) :: index, stride, upper
1046  sll_real64, dimension(:), intent(inout) :: data_out
1047 
1048  data_out(index) = data_out(index) + 0.5_f64*data_out(index - stride) &
1049  + 0.5_f64*data_out(modulo(index + stride - 1, upper) + 1);
1050  if (stride > 1) then
1051  call dehierarchical_stripe_recursive(index - stride/2, stride/2, upper, &
1052  data_out)
1053  call dehierarchical_stripe_recursive(index + stride/2, stride/2, upper, &
1054  data_out)
1055  end if
1056  end subroutine dehierarchical_stripe_recursive
1057 
1058 !!!! (De)hierarchization functions involving only parts of the dimensions
1059 
1060 ! Recursive worker function for dehierarchical_part
1061  recursive subroutine dehierarchical_part_d_dimension(interpolator,surplus,data_array,level,factor,index,dmax,dmin,dim,dorder)
1062  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
1063  sll_real64, intent(inout) :: surplus
1064  sll_real64, dimension(:), intent(in) :: data_array
1065  sll_int32, dimension(:), intent(in) :: level
1066  sll_real64, intent(in) :: factor
1067  sll_int32, intent(in) :: index
1068  sll_int32, intent(in) :: dmax, dmin, dim
1069  sll_int32, dimension(:), intent(in) :: dorder
1070  sll_int32 :: j, jj
1071 
1072  if (index .NE. -1) then
1073  if (dim > dmin - 1) then
1074  surplus = surplus + data_array(index)*factor
1075  end if
1076  do j = dim + 1, dmax
1077  jj = dorder(j)
1078  if (level(jj) > 0) then
1079  call dehierarchical_part_d_dimension(interpolator, surplus, data_array, &
1080  level, &
1081  -0.5_f64*factor, interpolator%hierarchy(index)%parent(2*jj - 1), &
1082  dmax, dmin, j, dorder)
1083  call dehierarchical_part_d_dimension(interpolator, surplus, data_array, &
1084  level, &
1085  -0.5_f64*factor, interpolator%hierarchy(index)%parent(2*jj), &
1086  dmax, dmin, j, dorder)
1087  end if
1088  end do
1089  end if
1090 
1091  end subroutine dehierarchical_part_d_dimension
1092 
1093 ! Recursive worker function for dehierarchical_part_order
1094 recursive subroutine dehierarchical_part_order_d_dimension(interpolator,surplus,data_array,start_level,level,k,factor,index,dmax,dmin,d,dorder)
1095  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
1096  sll_real64, intent(inout) :: surplus
1097  sll_real64, dimension(:), intent(in) :: data_array
1098  sll_int32, intent(in) :: start_level
1099  sll_int32, dimension(:), intent(in) :: level, k, dorder
1100  sll_real64, intent(in) :: factor
1101  sll_int32, intent(in) :: index
1102  sll_int32, intent(in) :: dmax, dmin, d
1103  sll_int32 :: j, father, jj
1104 
1105  if (d > dmin - 1) then
1106  surplus = surplus + data_array(index)*factor
1107  end if
1108  do j = d + 1, dmax
1109  jj = dorder(j)
1110  if (level(jj) > start_level) then
1111  father = max(interpolator%hierarchy(index)%parent(2*jj - 1), &
1112  interpolator%hierarchy(index)%parent(2*jj))
1113  if (father .NE. -1) then
1114  call dehierarchical_part_order_d_dimension( &
1115  interpolator, surplus, data_array, start_level, level, k, &
1116  interpolator%hs_weights(k(jj))*factor, father, dmax, dmin, j, dorder)
1117  end if
1118  end if
1119  end do
1120 
1121  end subroutine dehierarchical_part_order_d_dimension
1122 
1123 ! Dehierarchization from linear hierarchical surplus only applied to parts of the dimensions.
1124 ! Dimensions dehiearchized: dorder(dmin:dmax)
1125  subroutine dehierarchical_part(interpolator, data, dmax, dmin, dorder)
1126  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
1127  sll_real64, dimension(:), intent(inout) :: data
1128  sll_int32, intent(in) :: dmax, dmin
1129  sll_int32, dimension(:), intent(in) :: dorder
1130  !sll_real64, dimension(:), intent(out) :: data_out
1131  sll_int32 :: counter
1132  sll_real64 :: factor
1133 
1134  do counter = 1, interpolator%size_basis
1135  factor = -1.0_f64
1136  call dehierarchical_part_d_dimension &
1137  (interpolator, data(counter), &
1138  data, interpolator%hierarchy(counter)%level, &
1139  factor, counter, dmax, dmin, dmin - 1, dorder)
1140  end do
1141  end subroutine dehierarchical_part
1142 
1143 ! Dehierarchization from higher order hierarchical surplus only applied to parts of the dimensions.
1144 ! Dimensions dehiearchized: dorder(dmin:dmax)
1145  subroutine dehierarchical_part_order(interpolator, data, dmax, dmin, dorder, order)
1146  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
1147  sll_real64, dimension(:), intent(inout) :: data
1148  sll_int32, intent(in) :: dmax, dmin, order
1149  sll_int32, dimension(:), intent(in) :: dorder
1150  !sll_real64, dimension(:), intent(out) :: data_out
1151  sll_int32 :: counter, j, start_level, order_level_size
1152  sll_int32, dimension(:), allocatable :: k
1153  sll_real64 :: factor
1154 
1155  sll_allocate(k(interpolator%dim), j);
1156  start_level = order - 1
1157  order_level_size = 2**start_level
1158 
1159  do counter = 1, interpolator%size_basis
1160  factor = -1.0_f64
1161 
1162  do j = 1, interpolator%dim
1163  k(j) = modulo(interpolator%hierarchy(counter)%function_type(j) - &
1164  interpolator%hs_weights_index(order), &
1165  order_level_size) + interpolator%hs_weights_index(order);
1166  end do
1167  call dehierarchical_part_order_d_dimension &
1168  (interpolator, data(counter), &
1169  data, start_level, interpolator%hierarchy(counter)%level, k, &
1170  factor, counter, dmax, dmin, dmin - 1, dorder)
1171  end do
1172 
1173  end subroutine dehierarchical_part_order
1174 
1175 ! Hierarchization to linear hierarchical surplus only applied to parts of the dimensions.
1176 ! Dimensions hiearchized: dorder(dmin:dmax)
1177  subroutine hierarchical_part(interpolator, data, dmax, dmin, dorder)
1178  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
1179  sll_real64, dimension(:), intent(inout) :: data
1180  sll_int32, intent(in) :: dmax, dmin
1181  sll_int32, dimension(:), intent(in) :: dorder
1182  !sll_real64, dimension(:), intent(out) :: data_out
1183  sll_int32 :: counter
1184  sll_real64 :: factor
1185 
1186  do counter = interpolator%size_basis, 1, -1
1187  factor = 1.0_f64
1188  call dehierarchical_part_d_dimension &
1189  (interpolator, data(counter), &
1190  data, interpolator%hierarchy(counter)%level, &
1191  factor, counter, dmax, dmin, dmin - 1, dorder)
1192  end do
1193 
1194  end subroutine hierarchical_part
1195 
1196 ! Hierarchization to higher order hierarchical surplus only applied to parts of the dimensions.
1197 ! Dimensions hiearchized: dorder(dmin:dmax)
1198  subroutine hierarchical_part_order(interpolator, data, dmax, dmin, dorder, order)
1199  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
1200  sll_real64, dimension(:), intent(inout) :: data
1201  sll_int32, intent(in) :: dmax, dmin, order
1202  sll_int32, dimension(:), intent(in) :: dorder
1203  !sll_real64, dimension(:), intent(out) :: data_out
1204  sll_int32 :: counter, j
1205  sll_int32 :: start_level, order_level_size
1206  sll_int32, dimension(:), allocatable :: k
1207  sll_real64 :: factor
1208 
1209  sll_allocate(k(interpolator%dim), j);
1210  start_level = order - 1
1211  order_level_size = 2**start_level
1212 
1213  do counter = interpolator%size_basis, 1, -1
1214  do j = 1, interpolator%dim
1215  k(j) = modulo(interpolator%hierarchy(counter)%function_type(j) - &
1216  interpolator%hs_weights_index(order), &
1217  order_level_size) + interpolator%hs_weights_index(order);
1218  end do
1219  factor = 1.0_f64
1220  call dehierarchical_part_order_d_dimension &
1221  (interpolator, data(counter), &
1222  data, start_level, &
1223  interpolator%hierarchy(counter)%level, k, factor, counter, dmax, dmin, &
1224  dmin - 1, dorder)
1225  end do
1226  end subroutine hierarchical_part_order
1227 
1228 !!!! End general sparse grid helper functions !!!!
1229 !------------------------------------------------------------------------------!
1230 
1231 !------------------------------------------------------------------------------!
1232 !!!! Trapezoidal integrator on sparse grid
1233  subroutine integrate_trapezoidal(interpolator, data_in, val)
1234  class(sll_t_sparse_grid_interpolator), intent(in) :: interpolator
1235  sll_int32 :: i, j
1236  sll_real64, intent(inout) :: val
1237  sll_real64, dimension(:), intent(in) :: data_in
1238  sll_real64 :: phix1
1239 
1240  val = 0.0_f64
1241  do i = 0, interpolator%max_level
1242  if (interpolator%boundary == 0) then
1243  phix1 = 1.0_f64/(real((2**(i)), f64))
1244  else
1245  phix1 = 1.0_f64/(real(max(2**(i), 2), f64))
1246  end if
1247  do j = interpolator%level_mapping(i), interpolator%level_mapping(i + 1) - 1
1248  val = val + data_in(j)*phix1
1249  end do
1250  end do
1251  val = val*interpolator%volume
1252  end subroutine integrate_trapezoidal
1253 
1254 !Trapezoidal integrator on sparse grid working with the semi-hierarchical
1255 !surplus (hierarchical in (d-1) dimensions, nodal along one dimension). Note
1256 !this functions should usually not be used since it is only working if the
1257 !maximum number of levels along the nodal dimension is the maximum total levels
1258  subroutine integrate_trapezoidal2(interpolator, dorder, data_in, val)
1259  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
1260  sll_real64, intent(inout) :: val
1261  sll_real64, dimension(:), intent(inout) :: data_in
1262  sll_int32, dimension(:), intent(in) :: dorder
1263 
1264  call hierarchical_part(interpolator, data_in, interpolator%dim, 2, dorder)
1265  val = sum(data_in)*interpolator%volume/(real((2**(interpolator%levels(dorder(1)))), f64));
1266  call dehierarchical_part(interpolator, data_in, interpolator%dim, 2, dorder);
1267  end subroutine integrate_trapezoidal2
1268 
1269 !------------------------------------------------------------------------------!
1270 !!!!!! Various SGFFT helper functions. Need a clean up. !!!!
1271 
1272  subroutine extract_real_to_comp(sparsegrid, dim, max_level, index, data_in, data_out)
1273  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
1274  sll_int32, intent(in) :: dim, max_level, index
1275  sll_int32 :: n_points, index_running
1276  sll_real64, dimension(:), intent(in) :: data_in
1277  sll_comp64, dimension(:), intent(out) :: data_out
1278 
1279  n_points = 2**(max_level)
1280  data_out(1) = cmplx(data_in(index), 0.0_f64, kind=f64)
1281  if (max_level > 0) then
1282  index_running = sparsegrid%hierarchy(index)%children(dim*2)
1283 
1284  data_out(n_points/2 + 1) = cmplx(data_in(index_running), 0.0_f64, kind=f64)
1285  if (max_level > 1) then
1286 
1287  call extract_recursive_real_to_comp(sparsegrid, n_points/2 + 1, n_points/4, index_running, &
1288  dim, data_in, data_out)
1289  end if
1290  end if
1291 
1292  end subroutine extract_real_to_comp
1293 
1294  recursive subroutine extract_recursive_real_to_comp(sparsegrid, index_stripe, stride, index_sg, dim, data_in, data_out)
1295  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
1296  sll_int32, intent(in) :: index_stripe, stride, index_sg, dim
1297  sll_real64, dimension(:), intent(in) :: data_in
1298  sll_comp64, dimension(:), intent(inout) :: data_out
1299 
1300  data_out(index_stripe - stride) = &
1301  cmplx(data_in(sparsegrid%hierarchy(index_sg)%children(dim*2 - 1)), 0.0_f64, kind=f64)
1302  data_out(index_stripe + stride) = &
1303  cmplx(data_in(sparsegrid%hierarchy(index_sg)%children(dim*2)), 0.0_f64, kind=f64)
1304  if (stride > 1) then
1305  call extract_recursive_real_to_comp(sparsegrid, index_stripe - stride, stride/2, &
1306  sparsegrid%hierarchy(index_sg)%children(dim*2 - 1), dim, &
1307  data_in, data_out)
1308  call extract_recursive_real_to_comp(sparsegrid, index_stripe + stride, stride/2, &
1309  sparsegrid%hierarchy(index_sg)%children(dim*2), dim, &
1310  data_in, data_out)
1311  end if
1312  end subroutine extract_recursive_real_to_comp
1313 
1314  subroutine extract_comp(sparsegrid, dim, max_level, index, data_in, data_out)
1315  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
1316  sll_int32, intent(in) :: dim, max_level, index
1317  sll_int32 :: n_points, index_running
1318  sll_comp64, dimension(:), intent(in) :: data_in
1319  sll_comp64, dimension(:), intent(out) :: data_out
1320 
1321  n_points = 2**(max_level)
1322  data_out(1) = data_in(index)
1323  if (max_level > 0) then
1324  index_running = sparsegrid%hierarchy(index)%children(dim*2)
1325 
1326  data_out(n_points/2 + 1) = data_in(index_running)
1327  if (max_level > 1) then
1328 
1329  call extract_recursive_comp(sparsegrid, n_points/2 + 1, n_points/4, index_running, &
1330  dim, data_in, data_out)
1331  end if
1332  end if
1333 
1334  end subroutine extract_comp
1335 
1336  recursive subroutine extract_recursive_comp(sparsegrid, index_stripe, stride, index_sg, dim, data_in, data_out)
1337  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
1338  sll_int32, intent(in) :: index_stripe, stride, index_sg, dim
1339  sll_comp64, dimension(:), intent(in) :: data_in
1340  sll_comp64, dimension(:), intent(inout) :: data_out
1341 
1342  data_out(index_stripe - stride) = &
1343  data_in(sparsegrid%hierarchy(index_sg)%children(dim*2 - 1))
1344  data_out(index_stripe + stride) = &
1345  data_in(sparsegrid%hierarchy(index_sg)%children(dim*2))
1346  if (stride > 1) then
1347  call extract_recursive_comp(sparsegrid, index_stripe - stride, stride/2, &
1348  sparsegrid%hierarchy(index_sg)%children(dim*2 - 1), dim, &
1349  data_in, data_out)
1350  call extract_recursive_comp(sparsegrid, index_stripe + stride, stride/2, &
1351  sparsegrid%hierarchy(index_sg)%children(dim*2), dim, &
1352  data_in, data_out)
1353  end if
1354  end subroutine extract_recursive_comp
1355 
1356  subroutine extract_fourier(sparsegrid, dim, max_level, index, data_in, data_out)
1357  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
1358  sll_int32, intent(in) :: dim, max_level, index
1359  sll_int32 :: n_points, index_running
1360  sll_comp64, dimension(:), intent(in) :: data_in
1361  sll_comp64, dimension(:), intent(out) :: data_out
1362 
1363  n_points = 2**(max_level)
1364  data_out(1) = data_in(index)
1365  if (max_level > 0) then
1366  index_running = sparsegrid%hierarchy(index)%children(dim*2)
1367  data_out(2) = data_in(index_running)
1368  if (max_level > 1) then
1369  call extract_recursive_fourier(sparsegrid, index_running, 0, &
1370  2, max_level, dim, data_in, data_out)
1371  end if
1372  end if
1373  end subroutine extract_fourier
1374 
1375  recursive subroutine extract_recursive_fourier(sparsegrid, index_sg, ind, level, max_level, dim, data_in, data_out)
1376  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
1377  sll_int32, intent(in) :: level, max_level, index_sg, dim, ind
1378  sll_comp64, dimension(:), intent(in) :: data_in
1379  sll_comp64, dimension(:), intent(inout) :: data_out
1380 
1381  data_out(2**(level - 1) + 1 + 2*ind) = &
1382  data_in(sparsegrid%hierarchy(index_sg)%children(dim*2 - 1))
1383  data_out(2**(level - 1) + 1 + 2*ind + 1) = &
1384  data_in(sparsegrid%hierarchy(index_sg)%children(dim*2))
1385  if (level < max_level) then
1386  call extract_recursive_fourier(sparsegrid, &
1387  sparsegrid%hierarchy(index_sg)%children(dim*2 - 1), 2*ind, &
1388  level + 1, max_level, dim, data_in, data_out)
1389  call extract_recursive_fourier(sparsegrid, &
1390  sparsegrid%hierarchy(index_sg)%children(dim*2), 2*ind + 1, &
1391  level + 1, max_level, dim, data_in, data_out)
1392  end if
1393  end subroutine extract_recursive_fourier
1394 
1395  subroutine insert_fourier(sparsegrid, dim, max_level, index, data_in, data_out)
1396  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
1397  sll_int32, intent(in) :: dim, max_level, index
1398  sll_int32 :: n_points, index_running
1399  sll_comp64, dimension(:), intent(in) :: data_in
1400  sll_comp64, dimension(:), intent(out) :: data_out
1401 
1402  n_points = 2**(max_level)
1403  data_out(index) = data_in(1)
1404  if (max_level > 0) then
1405  index_running = sparsegrid%hierarchy(index)%children(dim*2)
1406  data_out(index_running) = data_in(2)
1407  if (max_level > 1) then
1408  call insert_recursive_fourier(sparsegrid, index_running, 0, &
1409  2, max_level, dim, data_in, data_out)
1410  end if
1411  end if
1412  end subroutine insert_fourier
1413 
1414  recursive subroutine insert_recursive_fourier(sparsegrid, index_sg, ind, level, max_level, dim, data_in, data_out)
1415  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
1416  sll_int32, intent(in) :: level, max_level, index_sg, dim, ind
1417  sll_comp64, dimension(:), intent(in) :: data_in
1418  sll_comp64, dimension(:), intent(inout) :: data_out
1419 
1420  data_out(sparsegrid%hierarchy(index_sg)%children(dim*2 - 1)) = &
1421  data_in(2**(level - 1) + 1 + 2*ind)
1422  data_out(sparsegrid%hierarchy(index_sg)%children(dim*2)) = &
1423  data_in(2**(level - 1) + 1 + 2*ind + 1)
1424  if (level < max_level) then
1425  call insert_recursive_fourier(sparsegrid, &
1426  sparsegrid%hierarchy(index_sg)%children(dim*2 - 1), 2*ind, &
1427  level + 1, max_level, dim, data_in, data_out)
1428  call insert_recursive_fourier(sparsegrid, &
1429  sparsegrid%hierarchy(index_sg)%children(dim*2), 2*ind + 1, &
1430  level + 1, max_level, dim, data_in, data_out)
1431  end if
1432  end subroutine insert_recursive_fourier
1433 
1434 !PN DEFINED BUT NOT USED
1435 !subroutine sort_hiera_1d(max_level,data_in,data_out)
1436 ! sll_int32, intent(in) :: max_level
1437 ! sll_real64,dimension(:),intent(in) :: data_in
1438 ! sll_real64,dimension(:),intent(inout) :: data_out
1439 ! sll_int32 :: size,l,j,initial,stride,counter
1440 !
1441 ! size = 2**max_level;
1442 !
1443 ! counter = 1
1444 ! data_out(1) = data_in(1);
1445 ! if(max_level>1) then
1446 ! data_out(2) = data_in(size/2+1)
1447 ! initial = size/4
1448 ! stride = size/2
1449 ! do l=2,max_level
1450 ! do j=0,2**(l-1)
1451 ! data_out(counter) = data_in(initial+1+j*stride)
1452 ! counter = counter +1
1453 ! end do
1454 ! initial = initial/2
1455 ! stride = stride/2
1456 ! end do
1457 ! end if
1458 !
1459 !end subroutine sort_hiera_1d
1460 
1461  subroutine insert_comp_to_real(sparsegrid, dim, max_level, index, data_in, data_out)
1462  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
1463  sll_int32, intent(in) :: dim, max_level, index
1464  sll_int32 :: n_points, index_running
1465  sll_comp64, dimension(:), intent(in) :: data_in
1466  sll_real64, dimension(:), intent(out) :: data_out
1467 
1468  n_points = 2**(max_level)
1469  data_out(index) = real(data_in(1))
1470  if (max_level > 0) then
1471  index_running = sparsegrid%hierarchy(index)%children(dim*2)
1472 
1473  data_out(index_running) = real(data_in(n_points/2 + 1))
1474  if (max_level > 1) then
1475 
1476  call insert_recursive_comp_to_real(sparsegrid, n_points/2 + 1, n_points/4, index_running, &
1477  dim, data_in, data_out)
1478  end if
1479  end if
1480 
1481  end subroutine insert_comp_to_real
1482 
1483  recursive subroutine insert_recursive_comp_to_real(sparsegrid, index_stripe, stride, index_sg, dim, data_in, data_out)
1484  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
1485  sll_int32, intent(in) :: index_stripe, stride, index_sg, dim
1486  sll_comp64, dimension(:), intent(in) :: data_in
1487  sll_real64, dimension(:), intent(inout) :: data_out
1488 
1489  data_out(sparsegrid%hierarchy(index_sg)%children(dim*2 - 1)) = &
1490  real(data_in(index_stripe - stride))
1491  data_out(sparsegrid%hierarchy(index_sg)%children(dim*2)) = &
1492  real(data_in(index_stripe + stride))
1493  if (stride > 1) then
1494  call insert_recursive_comp_to_real(sparsegrid, index_stripe - stride, stride/2, &
1495  sparsegrid%hierarchy(index_sg)%children(dim*2 - 1), dim, &
1496  data_in, data_out)
1497  call insert_recursive_comp_to_real(sparsegrid, index_stripe + stride, stride/2, &
1498  sparsegrid%hierarchy(index_sg)%children(dim*2), dim, &
1499  data_in, data_out)
1500  end if
1501  end subroutine insert_recursive_comp_to_real
1502 
1503  subroutine insert_comp(sparsegrid, dim, max_level, index, data_in, data_out)
1504  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
1505  sll_int32, intent(in) :: dim, max_level, index
1506  sll_int32 :: n_points, index_running
1507  sll_comp64, dimension(:), intent(in) :: data_in
1508  sll_comp64, dimension(:), intent(out) :: data_out
1509 
1510  n_points = 2**(max_level)
1511  data_out(index) = data_in(1)
1512  if (max_level > 0) then
1513  index_running = sparsegrid%hierarchy(index)%children(dim*2)
1514 
1515  data_out(index_running) = data_in(n_points/2 + 1)
1516  if (max_level > 1) then
1517 
1518  call insert_recursive_comp(sparsegrid, n_points/2 + 1, n_points/4, index_running, &
1519  dim, data_in, data_out)
1520  end if
1521  end if
1522 
1523  end subroutine insert_comp
1524 
1525  recursive subroutine insert_recursive_comp(sparsegrid, index_stripe, stride, index_sg, dim, data_in, data_out)
1526  class(sll_t_sparse_grid_interpolator), intent(in) :: sparsegrid
1527  sll_int32, intent(in) :: index_stripe, stride, index_sg, dim
1528  sll_comp64, dimension(:), intent(in) :: data_in
1529  sll_comp64, dimension(:), intent(inout) :: data_out
1530 
1531  data_out(sparsegrid%hierarchy(index_sg)%children(dim*2 - 1)) = &
1532  data_in(index_stripe - stride)
1533  data_out(sparsegrid%hierarchy(index_sg)%children(dim*2)) = &
1534  data_in(index_stripe + stride)
1535  if (stride > 1) then
1536  call insert_recursive_comp(sparsegrid, index_stripe - stride, stride/2, &
1537  sparsegrid%hierarchy(index_sg)%children(dim*2 - 1), dim, &
1538  data_in, data_out)
1539  call insert_recursive_comp(sparsegrid, index_stripe + stride, stride/2, &
1540  sparsegrid%hierarchy(index_sg)%children(dim*2), dim, &
1541  data_in, data_out)
1542  end if
1543  end subroutine insert_recursive_comp
1544 
1545  subroutine hira(data, max_level)
1546  sll_comp64, dimension(:), intent(inout) :: data
1547  sll_int32, intent(in) :: max_level
1548  sll_int32 ::l, i
1549 
1550  do l = max_level, 1, -1
1551  do i = 2**l, 2**(l - 1) + 1, -1
1552  data(2**l + 1 - i) = data(2**l + 1 - i) + data(i)
1553  end do
1554  end do
1555 
1556  end subroutine hira
1557 
1558  subroutine dehi(data, max_level)
1559  sll_comp64, dimension(:), intent(inout) :: data
1560  sll_int32, intent(in) :: max_level
1561  sll_int32 ::l, i
1562 
1563  do l = 1, max_level
1564  do i = 2**(l - 1) + 1, 2**l
1565  data(2**l + 1 - i) = data(2**l + 1 - i) - data(i)
1566  end do
1567  end do
1568 
1569  end subroutine dehi
1570 
1571  subroutine fft_on_stripe(interpolator, level)
1572  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
1573  sll_int32, intent(in) ::level
1574  sll_int32 :: i, no_points
1575 
1576  no_points = 2**level
1577 
1578  call sll_s_fft_exec_c2c_1d(interpolator%fft_object(level + 1)%fw, &
1579  interpolator%fft_object(level + 1)%in, interpolator%fft_object(level + 1)%out)
1580  ! Scale the Fourier coefficients
1581  do i = 1, no_points
1582  interpolator%fft_object(level + 1)%out(i) = &
1583  interpolator%fft_object(level + 1)%out(i)/cmplx(no_points, 0.0_f64, f64)
1584  end do
1585 
1586  end subroutine fft_on_stripe
1587 
1588  subroutine ifft_on_stripe(interpolator, level)
1589  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
1590  sll_int32, intent(in) :: level
1591 
1592  call sll_s_fft_exec_c2c_1d(interpolator%fft_object(level + 1)%bw, &
1593  interpolator%fft_object(level + 1)%out, interpolator%fft_object(level + 1)%in)
1594 
1595  end subroutine ifft_on_stripe
1596 
1597  subroutine fft_initialize(fft_object, levels)
1598  class(fft_hierarchical), dimension(:), pointer :: fft_object
1599  sll_int32, intent(in) :: levels
1600  sll_int32 :: l, size
1601 
1602  size = 1
1603  do l = 1, levels + 1
1604 
1605  fft_object(l)%sz_in = int(size, c_size_t)
1606  fft_object(l)%sz_out = int(size, c_size_t)
1607  fft_object(l)%out => sll_f_fft_allocate_aligned_complex(size)
1608  fft_object(l)%in => sll_f_fft_allocate_aligned_complex(size)
1609  call sll_s_fft_init_c2c_1d(fft_object(l)%fw, size, fft_object(l)%in, &
1610  fft_object(l)%out, sll_p_fft_backward, normalized=.false., &
1611  aligned=.true., optimization=sll_p_fft_measure)
1612 
1613  call sll_s_fft_init_c2c_1d(fft_object(l)%bw, size, fft_object(l)%out, &
1614  fft_object(l)%in, sll_p_fft_forward, normalized=.false., &
1615  aligned=.true., optimization=sll_p_fft_measure)
1616 
1617  size = size*2
1618  end do
1619 
1620  end subroutine fft_initialize
1621 
1622  subroutine fft_finalize(fft_object, levels)
1623  class(fft_hierarchical), dimension(:), pointer :: fft_object
1624  sll_int32, intent(in) :: levels
1625  sll_int32 :: l
1626 
1627  do l = 1, levels + 1
1628  call sll_s_fft_free(fft_object(l)%fw)
1629  call sll_s_fft_free(fft_object(l)%bw)
1630  end do
1631 
1632  end subroutine fft_finalize
1633 
1635  subroutine free_sparse_grid(interpolator)
1636  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
1637 
1638  call fft_finalize(interpolator%fft_object, interpolator%max_level)
1639  deallocate (interpolator%hierarchy)
1640  deallocate (interpolator%stripe)
1641  deallocate (interpolator%stripe_out)
1642  deallocate (interpolator%hs_weights)
1643  deallocate (interpolator%hs_weights_index)
1644  deallocate (interpolator%interp)
1645  deallocate (interpolator%level_mapping)
1646  deallocate (interpolator%eta_min)
1647  deallocate (interpolator%eta_max)
1648  deallocate (interpolator%length)
1649 
1650  end subroutine free_sparse_grid
1651 
1653  subroutine tohierarchical1d(interpolator, dim, max_level, index, data_in, data_out)
1654  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
1655  sll_real64, dimension(:), intent(in) :: data_in
1656  sll_comp64, dimension(:), intent(out) :: data_out
1657  sll_int32, intent(in) :: dim, max_level, index
1658 
1659  call extract_real_to_comp(interpolator, dim, max_level, &
1660  index, data_in, interpolator%fft_object(max_level + 1)%in)
1661  call fft_on_stripe(interpolator, max_level)
1662 
1663  call fft_to_centered(2**max_level, interpolator%fft_object(max_level + 1)%out, &
1664  interpolator%fft_object(max_level + 1)%in)
1665 
1666  call hira(interpolator%fft_object(max_level + 1)%in, max_level)
1667 
1668  call insert_fourier(interpolator, dim, max_level, &
1669  index, interpolator%fft_object(max_level + 1)%in, data_out)
1670  end subroutine tohierarchical1d
1671 
1673  subroutine tohierarchical1d_comp(interpolator, dim, max_level, index, data)
1674  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
1675  sll_comp64, dimension(:), intent(inout) :: data
1676  sll_int32, intent(in) :: dim, max_level, index
1677 
1678  call extract_comp(interpolator, dim, max_level, &
1679  index, data, interpolator%fft_object(max_level + 1)%in)
1680  call fft_on_stripe(interpolator, max_level)
1681 
1682  call fft_to_centered(2**max_level, interpolator%fft_object(max_level + 1)%out, &
1683  interpolator%fft_object(max_level + 1)%in)
1684  call hira(interpolator%fft_object(max_level + 1)%in, max_level)
1685  call insert_fourier(interpolator, dim, max_level, &
1686  index, interpolator%fft_object(max_level + 1)%in, data)
1687  end subroutine tohierarchical1d_comp
1688 
1690  subroutine tohira1d(interpolator, dim, max_level, index, data)
1691  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
1692  sll_comp64, dimension(:), intent(inout) :: data
1693  sll_int32, intent(in) :: dim, max_level, index
1694 
1695  call extract_fourier(interpolator, dim, max_level, index, &
1696  data, interpolator%fft_object(max_level + 1)%in)
1697  call hira(interpolator%fft_object(max_level + 1)%in, max_level)
1698  call insert_fourier(interpolator, dim, max_level, index, &
1699  interpolator%fft_object(max_level + 1)%in, data)
1700 
1701  end subroutine tohira1d
1702 
1704  subroutine tonodal1d(interpolator, dim, max_level, index, data_in, data_out)
1705  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
1706  sll_comp64, dimension(:), intent(in) :: data_in
1707  sll_real64, dimension(:), intent(out) :: data_out
1708  sll_int32, intent(in) :: dim, max_level, index
1709 
1710  call extract_fourier(interpolator, dim, &
1711  max_level, index, data_in, interpolator%fft_object(max_level + 1)%in)
1712  call dehi(interpolator%fft_object(max_level + 1)%in, max_level)
1713  call fft_to_inorder(2**max_level, interpolator%fft_object(max_level + 1)%in, &
1714  interpolator%fft_object(max_level + 1)%out)
1715  !call dehi(interpolator%fft_object(max_level+1)%out,max_level)
1716  call ifft_on_stripe(interpolator, max_level)
1717  call insert_comp_to_real(interpolator, dim, max_level, index, &
1718  interpolator%fft_object(max_level + 1)%in, data_out)
1719 
1720  end subroutine tonodal1d
1721 
1723  subroutine tonodal1d_comp(interpolator, dim, max_level, index, data)
1724  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
1725  sll_comp64, dimension(:), intent(inout) :: data
1726  sll_int32, intent(in) :: dim, max_level, index
1727 
1728  call extract_fourier(interpolator, dim, &
1729  max_level, index, data, interpolator%fft_object(max_level + 1)%in)
1730  call dehi(interpolator%fft_object(max_level + 1)%in, max_level)
1731  !print*, 'dehi', dim
1732  !print*, interpolator%fft_object(max_level+1)%in
1733  call fft_to_inorder(2**max_level, interpolator%fft_object(max_level + 1)%in, &
1734  interpolator%fft_object(max_level + 1)%out)
1735  !call dehi(interpolator%fft_object(max_level+1)%out,max_level)
1736  call ifft_on_stripe(interpolator, max_level)
1737  call insert_comp(interpolator, dim, max_level, index, &
1738  interpolator%fft_object(max_level + 1)%in, data)
1739 
1740  end subroutine tonodal1d_comp
1741 
1742  subroutine todehi1d(interpolator, dim, max_level, index, data_array)
1743  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
1744  sll_comp64, dimension(:), intent(inout) :: data_array
1745  sll_int32, intent(in) :: dim, max_level, index
1746 
1747  call extract_fourier(interpolator, dim, max_level, index, &
1748  data_array, interpolator%fft_object(max_level + 1)%out)
1749  call dehi(interpolator%fft_object(max_level + 1)%out, max_level)
1750  call insert_fourier(interpolator, dim, max_level, index, &
1751  interpolator%fft_object(max_level + 1)%out, data_array)
1752 
1753  end subroutine todehi1d
1754 
1755  subroutine fft_to_centered(length, data_in, data_out)
1756  sll_int32, intent(in) :: length
1757  sll_comp64, dimension(:), intent(in) :: data_in
1758  sll_comp64, dimension(:), intent(out) :: data_out
1759  sll_int32 :: k
1760 
1761  data_out(1) = data_in(1)
1762  data_out(length) = data_in(length/2 + 1)
1763 
1764  do k = 1, length/2 - 1
1765  data_out(2*k) = data_in(k + 1)
1766  data_out(2*k + 1) = data_in(length - k + 1)
1767  end do
1768 
1769  end subroutine fft_to_centered
1770 
1771  subroutine fft_to_inorder(length, data_in, data_out)
1772  sll_int32, intent(in) :: length
1773  sll_comp64, dimension(:), intent(in) :: data_in
1774  sll_comp64, dimension(:), intent(out) :: data_out
1775  sll_int32 :: k
1776 
1777  data_out(1) = data_in(1)
1778  data_out(length/2 + 1) = data_in(length)
1779 
1780  do k = 1, length/2 - 1
1781  data_out(k + 1) = data_in(2*k)
1782  data_out(length - k + 1) = data_in(2*k + 1)
1783  end do
1784 
1785  end subroutine fft_to_inorder
1786 
1787  subroutine displace_fourier_coeffs(d_scale, size, data)!displacement)
1788  sll_comp64, dimension(:), intent(inout) :: data
1789  !sll_real64, intent(in) :: displacement
1790  sll_real64, intent(in) :: d_scale
1791  sll_int32, intent(in) :: size
1792  sll_int32 :: k
1793  sll_real64 :: sin_val, cos_val
1794 
1795  do k = 1, size/2
1796  sin_val = sin(-real(k, f64)*d_scale)
1797  cos_val = cos(-real(k, f64)*d_scale)
1798  data(2*k) = cmplx(real(data(2*k))*cos_val - aimag(data(2*k))*sin_val, &
1799  real(data(2*k))*sin_val + aimag(data(2*k))*cos_val, kind=f64)
1800  if (k < size/2) then
1801  sin_val = -sin_val!sin(k*d_scale)
1802  !cos_val = cos(k*d_scale)
1803  data(2*k + 1) = cmplx(real(data(2*k + 1))*cos_val - aimag(data(2*k + 1))*sin_val, &
1804  real(data(2*k + 1))*sin_val + aimag(data(2*k + 1))*cos_val, kind=f64)
1805  end if
1806  end do
1807 
1808  end subroutine displace_fourier_coeffs
1809 
1810  subroutine displace1d(interpolator, dim, max_level, index, displacement, data)
1811  class(sll_t_sparse_grid_interpolator), intent(inout) :: interpolator
1812  sll_comp64, dimension(:), intent(inout) :: data
1813  sll_int32, intent(in) :: dim, max_level, index
1814  sll_real64, intent(in) :: displacement
1815 
1816  call extract_fourier(interpolator, dim, max_level, &
1817  index, data, interpolator%fft_object(max_level + 1)%in)
1818  call displace_fourier_coeffs(displacement, 2**max_level, &
1819  interpolator%fft_object(max_level + 1)%in)
1820  call insert_fourier(interpolator, dim, max_level, &
1821  index, interpolator%fft_object(max_level + 1)%in, data)
1822  end subroutine displace1d
1823 
1824 !!!! End SGFFT helper functions !!!!
1825 !------------------------------------------------------------------------------!
1826 
1827 end module sll_m_sparse_grid_interpolator
sll_m_fft
FFT interface for FFTW.
Definition: sll_m_fft_fftw3.F90:67
sll_m_fft::sll_s_fft_free
subroutine, public sll_s_fft_free(plan)
Delete a plan.
Definition: sll_m_fft_fftw3.F90:1007
sll_m_fft::sll_p_fft_forward
integer, parameter, public sll_p_fft_forward
Flag to specify forward FFT (negative sign) [value -1].
Definition: sll_m_fft_fftw3.F90:137
sll_m_fft::sll_s_fft_exec_c2c_1d
subroutine, public sll_s_fft_exec_c2c_1d(plan, array_in, array_out)
Compute fast Fourier transform in complex to complex mode.
Definition: sll_m_fft_fftw3.F90:414
sll_m_fft::sll_s_fft_init_c2c_1d
subroutine, public sll_s_fft_init_c2c_1d(plan, nx, array_in, array_out, direction, normalized, aligned, optimization)
Create new 1d complex to complex plan.
Definition: sll_m_fft_fftw3.F90:362
sll_m_fft::sll_p_fft_backward
integer, parameter, public sll_p_fft_backward
Flag to specify backward FFT (positive sign) [value 1].
Definition: sll_m_fft_fftw3.F90:138
sll_m_fft::sll_p_fft_measure
integer(c_int), parameter, public sll_p_fft_measure
FFTW planning-rigor flag FFTW_MEASURE (optimized plan) NOTE: planner overwrites the input array durin...
Definition: sll_m_fft_fftw3.F90:141
sll_m_fft::sll_f_fft_allocate_aligned_complex
complex(c_double_complex) function, dimension(:), pointer, public sll_f_fft_allocate_aligned_complex(n)
Function to allocate an aligned complex array.
Definition: sll_m_fft_fftw3.F90:178
sll_m_interpolators_1d_base
Module for 1D interpolation and reconstruction.
Definition: sll_m_interpolators_1d_base.F90:8
sll_m_periodic_interpolator_1d
Interpolator with periodic boundary conditions.
Definition: sll_m_periodic_interpolator_1d.F90:9
sll_m_sparse_grid_interpolator
Dimension-independent functions for sparse grid with polynomial basis functions.
Definition: sll_m_sparse_grid_interpolator.F90:6
sll_m_sparse_grid_interpolator::fft_to_inorder
subroutine fft_to_inorder(length, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:1772
sll_m_sparse_grid_interpolator::dehierarchical_part
subroutine dehierarchical_part(interpolator, data, dmax, dmin, dorder)
Definition: sll_m_sparse_grid_interpolator.F90:1126
sll_m_sparse_grid_interpolator::fft_on_stripe
subroutine fft_on_stripe(interpolator, level)
Definition: sll_m_sparse_grid_interpolator.F90:1572
sll_m_sparse_grid_interpolator::dehierarchical_stripe_order_recursive
recursive subroutine dehierarchical_stripe_order_recursive(index, stride, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:1031
sll_m_sparse_grid_interpolator::dehierarchical_part_d_dimension
recursive subroutine dehierarchical_part_d_dimension(interpolator, surplus, data_array, level, factor, index, dmax, dmin, dim, dorder)
Definition: sll_m_sparse_grid_interpolator.F90:1062
sll_m_sparse_grid_interpolator::extract_recursive_real_to_comp
recursive subroutine extract_recursive_real_to_comp(sparsegrid, index_stripe, stride, index_sg, dim, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:1295
sll_m_sparse_grid_interpolator::insert_recursive
recursive subroutine insert_recursive(sparsegrid, index_stripe, stride, index_sg, dim, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:714
sll_m_sparse_grid_interpolator::tonodal1d_comp
subroutine tonodal1d_comp(interpolator, dim, max_level, index, data)
Definition: sll_m_sparse_grid_interpolator.F90:1724
sll_m_sparse_grid_interpolator::dehierarchical_order_d_dimension
recursive subroutine dehierarchical_order_d_dimension(interpolator, surplus, data_array, start_level, level, k, factor, index, d)
Definition: sll_m_sparse_grid_interpolator.F90:916
sll_m_sparse_grid_interpolator::interpolate_disp_recursive
recursive subroutine interpolate_disp_recursive(interpolator, no_dims, dim, node_index, displacement, data_in, data_out, hiera)
Definition: sll_m_sparse_grid_interpolator.F90:505
sll_m_sparse_grid_interpolator::dehi
subroutine dehi(data, max_level)
Definition: sll_m_sparse_grid_interpolator.F90:1559
sll_m_sparse_grid_interpolator::displace_on_stripe_periodic
subroutine displace_on_stripe_periodic(interpolator, displacement, dim, max_level)
Definition: sll_m_sparse_grid_interpolator.F90:460
sll_m_sparse_grid_interpolator::hierarchical_stripe
subroutine hierarchical_stripe(sparsegrid, data, max_level)
Definition: sll_m_sparse_grid_interpolator.F90:947
sll_m_sparse_grid_interpolator::free_sparse_grid
subroutine free_sparse_grid(interpolator)
Finalize sparse grid.
Definition: sll_m_sparse_grid_interpolator.F90:1636
sll_m_sparse_grid_interpolator::insert_recursive_fourier
recursive subroutine insert_recursive_fourier(sparsegrid, index_sg, ind, level, max_level, dim, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:1415
sll_m_sparse_grid_interpolator::displace_fourier_coeffs
subroutine displace_fourier_coeffs(d_scale, size, data)
Definition: sll_m_sparse_grid_interpolator.F90:1788
sll_m_sparse_grid_interpolator::ifft_on_stripe
subroutine ifft_on_stripe(interpolator, level)
Definition: sll_m_sparse_grid_interpolator.F90:1589
sll_m_sparse_grid_interpolator::fft_to_centered
subroutine fft_to_centered(length, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:1756
sll_m_sparse_grid_interpolator::hierarchical
subroutine hierarchical(interpolator, data)
Definition: sll_m_sparse_grid_interpolator.F90:790
sll_m_sparse_grid_interpolator::initialize_sg
subroutine initialize_sg(interpolator, levels, order, interpolation, interpolation_type, eta_min, eta_max)
Definition: sll_m_sparse_grid_interpolator.F90:146
sll_m_sparse_grid_interpolator::extract_fourier
subroutine extract_fourier(sparsegrid, dim, max_level, index, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:1357
sll_m_sparse_grid_interpolator::hierarchical_order
subroutine hierarchical_order(interpolator, data, order)
Definition: sll_m_sparse_grid_interpolator.F90:808
sll_m_sparse_grid_interpolator::integrate_trapezoidal2
subroutine integrate_trapezoidal2(interpolator, dorder, data_in, val)
Definition: sll_m_sparse_grid_interpolator.F90:1259
sll_m_sparse_grid_interpolator::insert_periodic
subroutine insert_periodic(sparsegrid, dim, max_level, index, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:692
sll_m_sparse_grid_interpolator::dehierarchical_order
subroutine dehierarchical_order(interpolator, data, order)
Definition: sll_m_sparse_grid_interpolator.F90:856
sll_m_sparse_grid_interpolator::insert_periodic_additive
subroutine insert_periodic_additive(sparsegrid, factor, dim, max_level, index, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:734
sll_m_sparse_grid_interpolator::hierarchical_part
subroutine hierarchical_part(interpolator, data, dmax, dmin, dorder)
Definition: sll_m_sparse_grid_interpolator.F90:1178
sll_m_sparse_grid_interpolator::tonodal1d
subroutine tonodal1d(interpolator, dim, max_level, index, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:1705
sll_m_sparse_grid_interpolator::insert_recursive_comp
recursive subroutine insert_recursive_comp(sparsegrid, index_stripe, stride, index_sg, dim, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:1526
sll_m_sparse_grid_interpolator::compute_hierarchical_surplus
subroutine compute_hierarchical_surplus(interpolator, data_array)
Definition: sll_m_sparse_grid_interpolator.F90:283
sll_m_sparse_grid_interpolator::extract_recursive
recursive subroutine extract_recursive(sparsegrid, index_stripe, stride, index_sg, dim, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:669
sll_m_sparse_grid_interpolator::compute_linear_hierarchical_surplus
subroutine compute_linear_hierarchical_surplus(interpolator, data_array)
Definition: sll_m_sparse_grid_interpolator.F90:295
sll_m_sparse_grid_interpolator::insert_comp
subroutine insert_comp(sparsegrid, dim, max_level, index, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:1504
sll_m_sparse_grid_interpolator::extract_real_to_comp
subroutine extract_real_to_comp(sparsegrid, dim, max_level, index, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:1273
sll_m_sparse_grid_interpolator::tohierarchical1d_comp
subroutine tohierarchical1d_comp(interpolator, dim, max_level, index, data)
Complex version of ToHierarchical1d_comp.
Definition: sll_m_sparse_grid_interpolator.F90:1674
sll_m_sparse_grid_interpolator::hira
subroutine hira(data, max_level)
Definition: sll_m_sparse_grid_interpolator.F90:1546
sll_m_sparse_grid_interpolator::compute_dehierarchical
subroutine compute_dehierarchical(interpolator, data_array)
Definition: sll_m_sparse_grid_interpolator.F90:303
sll_m_sparse_grid_interpolator::interpolate_disp_1d_periodic_self
subroutine interpolate_disp_1d_periodic_self(interpolator, displacement, dim, max_level, index, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:628
sll_m_sparse_grid_interpolator::fft_finalize
subroutine fft_finalize(fft_object, levels)
Definition: sll_m_sparse_grid_interpolator.F90:1623
sll_m_sparse_grid_interpolator::basis_function
subroutine basis_function(x, fx, type)
Definition: sll_m_sparse_grid_interpolator.F90:327
sll_m_sparse_grid_interpolator::displace1d
subroutine displace1d(interpolator, dim, max_level, index, displacement, data)
Definition: sll_m_sparse_grid_interpolator.F90:1811
sll_m_sparse_grid_interpolator::insert_additive_recursive
recursive subroutine insert_additive_recursive(sparsegrid, factor, index_stripe, stride, index_sg, dim, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:757
sll_m_sparse_grid_interpolator::tohira1d
subroutine tohira1d(interpolator, dim, max_level, index, data)
Definition: sll_m_sparse_grid_interpolator.F90:1691
sll_m_sparse_grid_interpolator::basis_function_derivative
subroutine basis_function_derivative(x, fx, type)
Definition: sll_m_sparse_grid_interpolator.F90:397
sll_m_sparse_grid_interpolator::hierarchical_part_order
subroutine hierarchical_part_order(interpolator, data, dmax, dmin, dorder, order)
Definition: sll_m_sparse_grid_interpolator.F90:1199
sll_m_sparse_grid_interpolator::interpolate_disp_1d_periodic_for_neighbor
subroutine interpolate_disp_1d_periodic_for_neighbor(interpolator, displacement, factor, dim, max_level, max_level_neighbor, index, index_neighbor, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:602
sll_m_sparse_grid_interpolator::dehierarchical_part_order_d_dimension
recursive subroutine dehierarchical_part_order_d_dimension(interpolator, surplus, data_array, start_level, level, k, factor, index, dmax, dmin, d, dorder)
Definition: sll_m_sparse_grid_interpolator.F90:1095
sll_m_sparse_grid_interpolator::dehierarchical
subroutine dehierarchical(interpolator, data)
Definition: sll_m_sparse_grid_interpolator.F90:838
sll_m_sparse_grid_interpolator::extract_periodic
subroutine extract_periodic(sparsegrid, dim, max_level, index, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:649
sll_m_sparse_grid_interpolator::insert_fourier
subroutine insert_fourier(sparsegrid, dim, max_level, index, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:1396
sll_m_sparse_grid_interpolator::fft_initialize
subroutine fft_initialize(fft_object, levels)
Definition: sll_m_sparse_grid_interpolator.F90:1598
sll_m_sparse_grid_interpolator::interpolate_disp_1d_periodic
subroutine interpolate_disp_1d_periodic(interpolator, displacement, dim, max_level, index, data_in, data_out, hiera)
Definition: sll_m_sparse_grid_interpolator.F90:479
sll_m_sparse_grid_interpolator::extract_comp
subroutine extract_comp(sparsegrid, dim, max_level, index, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:1315
sll_m_sparse_grid_interpolator::displace_on_stripe_periodic_for_neighbor
subroutine displace_on_stripe_periodic_for_neighbor(interpolator, displacement, dim, max_level, max_level_neighbor)
Definition: sll_m_sparse_grid_interpolator.F90:441
sll_m_sparse_grid_interpolator::insert_recursive_comp_to_real
recursive subroutine insert_recursive_comp_to_real(sparsegrid, index_stripe, stride, index_sg, dim, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:1484
sll_m_sparse_grid_interpolator::extract_recursive_comp
recursive subroutine extract_recursive_comp(sparsegrid, index_stripe, stride, index_sg, dim, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:1337
sll_m_sparse_grid_interpolator::integrate_trapezoidal
subroutine integrate_trapezoidal(interpolator, data_in, val)
Definition: sll_m_sparse_grid_interpolator.F90:1234
sll_m_sparse_grid_interpolator::dehierarchical_stripe
subroutine dehierarchical_stripe(sparsegrid, data, max_level)
Definition: sll_m_sparse_grid_interpolator.F90:991
sll_m_sparse_grid_interpolator::interpolate_disp
subroutine interpolate_disp(interpolator, dim, displacement, data_in, data_out, hiera)
Definition: sll_m_sparse_grid_interpolator.F90:539
sll_m_sparse_grid_interpolator::dehierarchical_d_dimension
recursive subroutine dehierarchical_d_dimension(interpolator, surplus, data_array, level, factor, index, d)
Definition: sll_m_sparse_grid_interpolator.F90:886
sll_m_sparse_grid_interpolator::insert_comp_to_real
subroutine insert_comp_to_real(sparsegrid, dim, max_level, index, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:1462
sll_m_sparse_grid_interpolator::dehierarchical_part_order
subroutine dehierarchical_part_order(interpolator, data, dmax, dmin, dorder, order)
Definition: sll_m_sparse_grid_interpolator.F90:1146
sll_m_sparse_grid_interpolator::tohierarchical1d
subroutine tohierarchical1d(interpolator, dim, max_level, index, data_in, data_out)
Compute Fourier coefficients on sparse grid along dimension dim.
Definition: sll_m_sparse_grid_interpolator.F90:1654
sll_m_sparse_grid_interpolator::dehierarchical_stripe_recursive
recursive subroutine dehierarchical_stripe_recursive(index, stride, upper, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:1045
sll_m_sparse_grid_interpolator::todehi1d
subroutine todehi1d(interpolator, dim, max_level, index, data_array)
Definition: sll_m_sparse_grid_interpolator.F90:1743
sll_m_sparse_grid_interpolator::extract_recursive_fourier
recursive subroutine extract_recursive_fourier(sparsegrid, index_sg, ind, level, max_level, dim, data_in, data_out)
Definition: sll_m_sparse_grid_interpolator.F90:1376
sll_m_fft::sll_t_fft
Type for FFT plan in SeLaLib.
Definition: sll_m_fft_fftw3.F90:123
sll_m_interpolators_1d_base::sll_c_interpolator_1d
Abstract class for 1D interpolation and reconstruction.
Definition: sll_m_interpolators_1d_base.F90:21
sll_m_periodic_interpolator_1d::sll_t_periodic_interpolator_1d
Periodic interpolator.
Definition: sll_m_periodic_interpolator_1d.F90:34
sll_m_sparse_grid_interpolator::fft_hierarchical
class to hold values for hierarchical fft computations
Definition: sll_m_sparse_grid_interpolator.F90:41
sll_m_sparse_grid_interpolator::interpolator_base_ptr
Definition: sll_m_sparse_grid_interpolator.F90:62
sll_m_sparse_grid_interpolator::sll_t_sparse_grid_interpolator
Class defining the sparse grid data structure.
Definition: sll_m_sparse_grid_interpolator.F90:67
sll_m_sparse_grid_interpolator::sparsegrid_node
Data type for sparse grid node.
Definition: sll_m_sparse_grid_interpolator.F90:51
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:39 for Semi-Lagrangian Library by   doxygen 1.9.1