Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_sparse_grid_2d.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_sg2d
33  procedure :: interpolate_from_interpolant_value => interpolate_value_sg
36  procedure :: fg_to_sg
37  procedure :: spfft
38  procedure :: ispfft
40  ! Filter functions that can be used for stabilization
41  procedure :: filter_highest
42  procedure :: filter
43  procedure :: linear_filter
44  ! From here these are helper functions
45  procedure :: displace
46  procedure :: set_hierarchy_info
48  procedure :: tohira
49  procedure :: todehi
50  procedure :: tohierarchical
51  procedure :: tonodal
52  end type sll_t_sparse_grid_2d
53 
54 contains
55 
56 !------------------------------------------------------------------------------!
57 !!!! SGFFT routines !!!!!
58 
60  subroutine interpolate_array_disp_sgfft(interpolator, dim, displacment_in, data_in, data_out)
61  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
62  sll_int32, intent(in) :: dim
63  sll_real64, intent(in) :: displacment_in
64  sll_comp64, dimension(:), intent(inout) :: data_in
65  sll_real64, dimension(:), intent(out) :: data_out
66 
67  sll_real64:: displacement
68 
69  displacement = displacment_in*2.0_f64*sll_p_pi/interpolator%length(dim)
70  call interpolator%displace(dim, displacement, data_in);
71  call interpolator%ISPFFT(data_in, data_out);
72  end subroutine interpolate_array_disp_sgfft
73 
74 !!!! End SGFFT routines !!!!!
75 !------------------------------------------------------------------------------!
76 
77 !------------------------------------------------------------------------------!
78 !!!!!!!!!! Initialization routines !!!!!!!!
79 
81  subroutine initialize_sg2d( &
82  interpolator, &
83  levels, &
84  order, &
85  interpolation, &
86  interpolation_type, &
87  eta_min, &
88  eta_max, &
89  boundary, &
90  modified)
91 
92  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
93  sll_real64, dimension(:), intent(in) :: eta_min
94  sll_real64, dimension(:), intent(in) :: eta_max
95  sll_int32, dimension(:), intent(in) :: levels
96  sll_int32, intent(in) :: order
97  sll_int32, intent(in) :: interpolation
98  sll_int32, intent(in) :: interpolation_type
99  sll_int32, intent(in) :: modified
100 
103  sll_int32, intent(in) :: boundary
104 
105  sll_int32 :: i, j, k1, k2, l1, l2, l, counter
106  sll_int32 :: ierr
107  sll_int32, dimension(:), allocatable :: novec, lvec, kvec
108  sll_int32 :: no1, no2
109 
110  interpolator%dim = 2;
111  interpolator%modified = modified;
112  interpolator%boundary = boundary;
113  interpolator%max_level = levels(1);
114  do l = 2, interpolator%dim
115  interpolator%max_level = max(levels(l), interpolator%max_level);
116  end do
117  if (interpolator%modified == 1) then
118  interpolator%max_level = interpolator%max_level + 1;
119  end if
120 
121  interpolator%size_basis = 0;
122  if (interpolator%boundary == 0) then
123  do l = 0, interpolator%max_level
124  do l1 = max(0, l - levels(2)), min(l, levels(1))
125  l2 = l - l1
126  interpolator%size_basis = &
127  interpolator%size_basis + &
128  max(2**(l1 - 1), 1)*max(2**(l2 - 1), 1);
129  end do
130  end do
131  else
132  do l = 0, interpolator%max_level
133  do l1 = max(0, l - levels(2)), min(l, levels(1))
134  l2 = l - l1
135  if (l1 == 0) then
136  no1 = 2;
137  else
138  no1 = 2**(l1 - 1);
139  end if
140  if (l2 == 0) then
141  no2 = 2;
142  else
143  no2 = 2**(l2 - 1);
144  end if
145  interpolator%size_basis = &
146  interpolator%size_basis + no1*no2;
147  end do
148  end do
149  end if
150 
151  call interpolator%initialize_sg(levels, order, interpolation, &
152  interpolation_type, eta_min, eta_max);
153  sll_allocate(lvec(interpolator%dim), ierr);
154  sll_allocate(kvec(interpolator%dim), ierr);
155  sll_allocate(novec(interpolator%dim), ierr);
156  sll_allocate(interpolator%index(0:interpolator%levels(1), 0:interpolator%levels(2)), ierr);
157  ! Set the hierarchy of the grid
158  counter = 1
159  do l = 0, interpolator%max_level
160  interpolator%level_mapping(l) = counter;
161  do l1 = max(0, l - levels(2)), min(l, levels(1))
162  novec(1) = max(2**(l1 - 1), 1)
163  if (interpolator%boundary == 1 .AND. l1 == 0) then
164  novec(1) = 2;
165  end if
166  lvec(1) = l1;
167  l2 = l - l1
168  novec(2) = max(2**(l2 - 1), 1)
169  if (interpolator%boundary == 1 .AND. l2 == 0) then
170  novec(2) = 2;
171  end if
172  lvec(2) = l2;
173  interpolator%index(l1, l2) = counter
174  do k1 = 0, novec(1) - 1
175  kvec(1) = k1
176  do k2 = 0, novec(2) - 1
177  kvec(2) = k2
178  interpolator%hierarchy(counter)%index_on_level(1) = k1;
179  interpolator%hierarchy(counter)%index_on_level(2) = k2;
180  do j = 1, interpolator%dim
181  if (interpolator%boundary == 0) then
182  call interpolator%set_hierarchy_info &
183  (counter, j, lvec, kvec, novec);
184  elseif (interpolator%boundary == 1) then
185  call interpolator%set_hierarchy_info_boundary &
186  (counter, j, lvec, kvec, novec);
187  end if
188  end do
189  counter = counter + 1;
190  end do
191  end do
192  end do
193  end do
194  interpolator%level_mapping(interpolator%max_level + 1) = counter;
195  ! Now rescale all the coordinates to the actual mesh size
196  do i = 1, interpolator%size_basis
197  do j = 1, interpolator%dim
198  interpolator%hierarchy(i)%coordinate(j) = &
199  interpolator%hierarchy(i)%coordinate(j)* &
200  interpolator%length(j) + &
201  interpolator%eta_min(j)
202  end do
203  end do
204 
205  end subroutine initialize_sg2d
206 
208  subroutine set_hierarchy_info(interpolator, counter, cdim, lvecin, kvecin, novecin)
209  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
210  sll_int32 :: ld
211  sll_int32 :: kd
212  sll_int32, intent(in) :: cdim
213  sll_int32, intent(in) :: counter ! counter for node
214  sll_int32, dimension(:), intent(in) :: lvecin
215  sll_int32, dimension(:), intent(in) :: kvecin
216  sll_int32, dimension(:), intent(in) :: novecin
217  sll_int32, dimension(:), allocatable :: lvec, kvec, novec
218  sll_int32 :: jj, stride
219 
220  sll_allocate(lvec(interpolator%dim), jj);
221  sll_allocate(kvec(interpolator%dim), jj);
222  sll_allocate(novec(interpolator%dim), jj);
223  do jj = 1, interpolator%dim
224  lvec(jj) = lvecin(jj);
225  kvec(jj) = kvecin(jj);
226  novec(jj) = novecin(jj);
227  end do
228  ld = lvec(cdim);
229  kd = kvec(cdim);
230  interpolator%hierarchy(counter)%level(cdim) = ld;
231  stride = cdim*2 - 1
232  if (ld == 0) then
233  interpolator%hierarchy(counter)%coordinate(cdim) = 0.0_f64
234  interpolator%hierarchy(counter)%parent(stride) = &
235  counter
236  interpolator%hierarchy(counter)%parent(stride + 1) = &
237  counter
238  interpolator%hierarchy(counter)%function_type(cdim) = 0
239  else
240  interpolator%hierarchy(counter)%coordinate(cdim) = &
241  1.0_f64/(2.0_f64**ld) + real(kd, f64)*1.0_f64/(2.0_f64**(ld - 1))
242 
243  lvec(cdim) = lvec(cdim) - 1;
244  novec(cdim) = max(novec(cdim)/2, 1);
245  ! This one is actually only a neighboring point not the direct parent.
246  kvec(cdim) = modulo((kd + 1)/2, max(2**(ld - 2), 1));
247  interpolator%hierarchy(counter)%parent( &
248  modulo(kd, 2) + stride) = &
249  interpolator%hierarchy( &
250  interpolator%index(lvec(1), lvec(2)) + &
251  kvec(1)*novec(2) + kvec(2))%parent(stride)
252  ! This is the actual parent.
253  kvec(cdim) = kd/2;
254  interpolator%hierarchy(counter)%parent( &
255  modulo(kd + 1, 2) + stride) = &
256  interpolator%index(lvec(1), lvec(2)) + &
257  kvec(1)*novec(2) + kvec(2)
258  ! Now tell my parent that I am his child.
259  interpolator%hierarchy( &
260  interpolator%hierarchy(counter)%parent( &
261  modulo(kd + 1, 2) + stride))%children(modulo(kd, 2) + stride) = counter
262  if (ld == 1) then
263  interpolator%hierarchy( &
264  interpolator%hierarchy(counter)%parent( &
265  modulo(kd + 1, 2) + stride))%children(modulo(kd + 1, 2) + stride) = counter
266  end if
267  if (interpolator%order == 1) then
268  interpolator%hierarchy(counter)% &
269  function_type(cdim) = -1
270  elseif (ld == 1 .OR. interpolator%order == 2) then
271  interpolator%hierarchy(counter)% &
272  function_type(cdim) = 1 + modulo(kd, 2)
273  elseif (ld == 2 .OR. interpolator%order == 3) then !(order==3) then!
274  interpolator%hierarchy(counter)% &
275  function_type(cdim) = 3 + modulo(kd, 4)
276  else
277  interpolator%hierarchy(counter)% &
278  function_type(cdim) = 7 + modulo(kd, 8)
279  end if
280  end if
281 
282  end subroutine set_hierarchy_info
283 
285  subroutine set_hierarchy_info_boundary(interpolator, counter, cdim, lvecin, kvecin, novecin)
286  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
287  sll_int32 :: ld
288  sll_int32 :: kd
289  sll_int32, intent(in) :: cdim
290  sll_int32, intent(in) :: counter
291  sll_int32, dimension(:), intent(in) :: lvecin, kvecin, novecin
292  sll_int32, dimension(:), allocatable :: lvec, kvec, novec
293  sll_int32 :: jj, stride
294 
295  sll_allocate(lvec(interpolator%dim), jj);
296  sll_allocate(kvec(interpolator%dim), jj);
297  sll_allocate(novec(interpolator%dim), jj);
298  do jj = 1, interpolator%dim
299  lvec(jj) = lvecin(jj);
300  kvec(jj) = kvecin(jj);
301  novec(jj) = novecin(jj);
302  end do
303  ld = lvec(cdim);
304  kd = kvec(cdim);
305  interpolator%hierarchy(counter)%level(cdim) = ld;
306  stride = cdim*2 - 1
307  if (ld == 0) then
308  interpolator%hierarchy(counter)%coordinate(cdim) = real(kd, f64);
309  interpolator%hierarchy(counter)%parent(stride) = &
310  counter
311  interpolator%hierarchy(counter)%parent(stride + 1) = &
312  counter
313  interpolator%hierarchy(counter)%function_type(cdim) = -1
314  else
315  interpolator%hierarchy(counter)%coordinate(cdim) = &
316  1.0_f64/(2.0_f64**ld) + real(kd, f64)*1.0_f64/(2.0_f64**(ld - 1))
317 
318  lvec(cdim) = lvec(cdim) - 1;
319  novec(cdim) = novec(cdim)/2;
320  if (ld == 1) then
321  kvec(cdim) = kd/2;
322  novec(cdim) = 2;
323  interpolator%hierarchy(counter)%parent(stride) = &
324  interpolator%index(lvec(1), lvec(2)) + &
325  kvec(1)*novec(2) + kvec(2)
326  if (cdim == 1) then
327  interpolator%hierarchy(counter)%parent( &
328  stride + 1) = interpolator%hierarchy(counter)%parent( &
329  stride) + novec(2)
330  else
331  interpolator%hierarchy(counter)%parent( &
332  stride + 1) = interpolator%hierarchy(counter)%parent( &
333  stride) + 1
334  end if
335  interpolator%hierarchy(interpolator%hierarchy(counter)%parent(stride))%children &
336  (stride) = counter
337  !interpolator%hierarchy(interpolator%hierarchy(counter)%parent(stride+1))%&
338  ! children(stride) = counter
339  else
340  ! This one is actually only a neighboring point not the direct parent.
341  kvec(cdim) = modulo((kd + 1)/2, max(2**(ld - 2), 1));
342  kvec(cdim) = (kd + 1)/2;
343  if (kvec(cdim) < novec(cdim)) then
344  interpolator%hierarchy(counter)%parent( &
345  modulo(kd, 2) + stride) = &
346  interpolator%hierarchy( &
347  interpolator%index(lvec(1), lvec(2)) + &
348  kvec(1)*novec(2) + kvec(2))%parent(stride)
349  else
350  kvec(cdim) = kvec(cdim) - 1;
351  interpolator%hierarchy(counter)%parent( &
352  modulo(kd, 2) + stride) = &
353  interpolator%hierarchy( &
354  interpolator%index(lvec(1), lvec(2)) + &
355  kvec(1)*novec(2) + kvec(2))%parent(stride + 1)
356  end if
357  ! This is the actual parent.
358  kvec(cdim) = kd/2;
359  interpolator%hierarchy(counter)%parent( &
360  modulo(kd + 1, 2) + stride) = &
361  interpolator%index(lvec(1), lvec(2)) + &
362  kvec(1)*novec(2) + kvec(2)
363  ! Now tell my parent that I am his child.
364  interpolator%hierarchy( &
365  interpolator%hierarchy(counter)%parent( &
366  modulo(kd + 1, 2) + stride))%children(modulo(kd, 2) + stride) = counter
367  end if
368 
369  if (interpolator%order == 1) then
370  interpolator%hierarchy(counter)% &
371  function_type(cdim) = -1
372  elseif (ld == 1 .OR. interpolator%order == 2) then
373  interpolator%hierarchy(counter)% &
374  function_type(cdim) = 1 + modulo(kd, 2)
375  elseif (ld == 2 .OR. interpolator%order == 3) then !(order==3) then!
376  interpolator%hierarchy(counter)% &
377  function_type(cdim) = 3 + modulo(kd, 4)
378  else
379  interpolator%hierarchy(counter)% &
380  function_type(cdim) = 7 + modulo(kd, 8)
381  end if
382  end if
383 
384  end subroutine set_hierarchy_info_boundary
385 
386 !!!! End initialization routines !!!!
387 !------------------------------------------------------------------------------!
388 
389 !------------------------------------------------------------------------------!
390 !!!!! Interpolation functions !!!!!!
391 
393  function interpolate_value_sg(interpolator, data, eta) result(val)
394  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
395  sll_real64, dimension(:), intent(in) :: data
396  sll_real64 :: val
397  sll_real64, dimension(:), intent(in) :: eta
398 
399  if (interpolator%boundary == 0) then
400  val = interpolate_from_hierarchical_surplus(interpolator, data, eta)
401  else
402  val = interpolate_from_hierarchical_surplus_boundary(interpolator, data, eta)
403  end if
404 
405  end function interpolate_value_sg
406 
407 ! helper function for interpolate_from_interpolant_value
408 
410  function interpolate_from_hierarchical_surplus(interpolator, data, eta) result(val)
411  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
412  sll_int32 :: j, l1, l2, level
413  sll_real64 :: val
414  sll_real64, dimension(:), intent(in) :: data
415  sll_real64, dimension(:), intent(in) :: eta
416  sll_real64, dimension(2) :: eta_norm
417  sll_real64, dimension(2) :: phi
418  sll_int32, dimension(2) :: no, l
419  sll_int32, dimension(:, :), allocatable :: ind
420  sll_real64 :: scale
421  sll_int32 :: index
422 
423  sll_allocate(ind(0:interpolator%max_level, 1:interpolator%dim), j)
424 
425  val = 0.0_f64
426  ind(0:1, 1:interpolator%dim) = 0
427 
428  do j = 1, interpolator%dim
429  eta_norm(j) = (eta(j) - interpolator%eta_min(j))/interpolator%length(j)
430  eta_norm(j) = modulo(eta_norm(j), 1.0_f64)
431 
432  scale = 0.5_f64
433  do level = 2, interpolator%levels(j)
434  ind(level, j) = ind(level - 1, j)*2
435  if (eta_norm(j) > scale*real(ind(level, j) + 1, f64)) then
436  ind(level, j) = ind(level, j) + 1
437  end if
438  scale = scale*0.5_f64
439  end do
440  end do
441 
442  do level = 0, interpolator%max_level
443  do l1 = max(0, level - interpolator%levels(2)), min(level, interpolator%levels(1))
444  l(1) = l1
445  no(1) = max(2**(l1 - 1), 1)
446  l2 = level - l1
447  l(2) = l2
448  no(2) = max(2**(l2 - 1), 1)
449  index = interpolator%index(l1, l2) + ind(l1, 1)*no(2) &
450  + ind(l2, 2)
451  do j = 1, interpolator%dim
452  call interpolator%basis_function(real(2**(max(l(j), 1)), f64)*eta_norm(j) &
453  - real(2*ind(l(j), j), f64) - 1.0_f64, phi(j), &
454  interpolator%hierarchy(index)%function_type(j))
455  end do
456  val = val + data(index) &
457  *phi(1)*phi(2)
458  end do
459  end do
460 
462 
463 ! interpolation from hierarchical surplus non-periodic
464 
466  function interpolate_from_hierarchical_surplus_boundary(interpolator, data, eta) result(val)
467  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
468  sll_int32 :: j, l1, l2, level
469  sll_real64 :: val
470  sll_real64, dimension(:), intent(in) :: data
471  sll_real64, dimension(:), intent(in) :: eta
472  sll_real64, dimension(2) :: eta_norm
473  sll_real64, dimension(2) :: phi
474  sll_int32, dimension(2) :: no, l
475  sll_int32, dimension(:, :), allocatable :: ind
476  sll_real64 :: scale, phi1a, phi2a
477  sll_int32 :: index
478 
479  sll_allocate(ind(0:interpolator%max_level, 1:interpolator%dim), j)
480 
481  val = 0.0_f64
482  ind(0:1, 1:interpolator%dim) = 0
483 
484  do j = 1, interpolator%dim
485  eta_norm(j) = (eta(j) - interpolator%eta_min(j))/interpolator%length(j)
486 
487  scale = 0.5_f64
488  do level = 2, interpolator%levels(j)
489  ind(level, j) = ind(level - 1, j)*2
490  if (eta_norm(j) > scale*real(ind(level, j) + 1, f64)) then
491  ind(level, j) = ind(level, j) + 1
492  end if
493  scale = scale*0.5_f64
494  end do
495  end do
496 
497  do level = 0, interpolator%max_level
498  do l1 = max(0, level - interpolator%levels(2)), min(level, interpolator%levels(1))
499  l(1) = l1
500  if (l1 == 0) then
501  no(1) = 2;
502  else
503  no(1) = 2**(l1 - 1);
504  end if
505  l2 = level - l1
506  l(2) = l2
507  if (l2 == 0) then
508  no(2) = 2;
509  else
510  no(2) = max(2**(l2 - 1), 1);
511  end if
512  if (l1 > 0) then
513  if (l2 > 0) then
514  index = interpolator%index(l1, l2) + ind(l1, 1)*no(2) &
515  + ind(l2, 2)
516 
517  do j = 1, interpolator%dim
518  call interpolator%basis_function(real(2**l(j), f64)*eta_norm(j) &
519  - real(2*ind(l(j), j), f64) - 1.0_f64, phi(j), &
520  interpolator%hierarchy(index)%function_type(j))
521  end do
522  val = val + data(index) &
523  *phi(1)*phi(2)
524  else ! l2=0
525  index = interpolator%index(l1, l2) + ind(l1, 1)*no(2)
526  call interpolator%basis_function(real(2**l(1), f64)*eta_norm(1) &
527  - real(2*ind(l(1), 1), f64) - 1.0_f64, phi(1), &
528  interpolator%hierarchy(index)%function_type(1))
529  call interpolator%basis_function(eta_norm(2), phi(2), -1)
530  val = val + data(index)*phi(1)*phi(2);
531  call interpolator%basis_function(eta_norm(2) - 1.0_f64, phi(2), -1)
532  val = val + data(index + 1)*phi(1)*phi(2);
533  end if
534  else !l1 = 0
535  if (l2 > 0) then
536  index = interpolator%index(l1, l2) + ind(l2, 2)
537  call interpolator%basis_function(real(2**l(2), f64)*eta_norm(2) &
538  - real(2*ind(l(2), 2), f64) - 1.0_f64, phi(2), &
539  interpolator%hierarchy(index)%function_type(2))
540  call interpolator%basis_function(eta_norm(1), phi(1), -1)
541  val = val + data(index)*phi(1)*phi(2);
542  call interpolator%basis_function(eta_norm(1) - 1.0_f64, phi(1), -1)
543  val = val + data(index + no(2))*phi(1)*phi(2);
544  else !l1=l2=0
545  index = interpolator%index(l1, l2);
546  call interpolator%basis_function(eta_norm(1), phi(1), -1)
547  call interpolator%basis_function(eta_norm(1) - 1.0_f64, phi1a, -1)
548  call interpolator%basis_function(eta_norm(2), phi(2), -1)
549  call interpolator%basis_function(eta_norm(2) - 1.0_f64, phi2a, -1)
550  val = val + data(index)*phi(1)*phi(2) + &
551  data(index + 1)*phi(1)*phi2a + data(index + 2)*phi1a*phi(2) + &
552  data(index + 3)*phi1a*phi2a
553  end if
554  end if
555  end do
556  end do
557 
559 
565  subroutine interpolate_const_disp(interpolator, dorder, displacement, data_in, data_out, hiera)
566  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
567  sll_int32, dimension(:), intent(in) :: dorder
568  sll_real64, dimension(:), intent(inout) :: data_in
569  sll_real64, dimension(:), intent(out) :: data_out
570  sll_real64, intent(in) ::displacement
571  logical, intent(in) :: hiera
572 
573  sll_int32 :: j, counter, i2, k2
574  sll_int32, dimension(2) :: ind_order, no
575  sll_int32, dimension(:, :), allocatable :: ind
576 
577  sll_allocate(ind(interpolator%max_level + 1, 2), i2);
578  ind_order(dorder(1)) = 0
579  no(dorder(1)) = 1
580  do i2 = 0, interpolator%levels(dorder(2))
581  ind_order(dorder(2)) = i2
582  no(dorder(2)) = max(2**(i2 - 1), 1);
583  ind(1, dorder(1)) = 0;
584  do k2 = 0, no(dorder(2)) - 1
585  ind(ind_order(dorder(2)) + 1, dorder(2)) = k2;
586  counter = interpolator%index( &
587  ind_order(1), ind_order(2)) + &
588  ind(ind_order(1) + 1, 1)*no(2) &
589  + ind(ind_order(2) + 1, 2)
590 
591  ! Evaluate along dorder(1)-stripe
592  call interpolator%interpolate_disp_1d_periodic &
593  (displacement, dorder(1), &
594  min(interpolator%levels(dorder(1)), interpolator%max_level &
595  - ind_order(dorder(2))), counter, data_in, data_out, hiera);
596  end do
597  end do
598 
599  if (hiera .EQV. .false.) then
600  ! Dehierarchization along dimension dorder(1) only
601  do j = interpolator%order, 2, -1
602  call interpolator%dehierarchical_part_order &
603  (data_out, &
604  interpolator%dim, 2, dorder, j)
605  end do
606 
607  call interpolator%dehierarchical_part(data_out, &
608  interpolator%dim, 2, dorder)
609  end if
610 
611  end subroutine interpolate_const_disp
612 
617 
618  subroutine interpolate_disp_nconst_in_1d(interpolator, displacement, dorder, data_in, data_out)
619  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
620  sll_real64, dimension(:), intent(inout) :: data_in
621  sll_real64, dimension(:), intent(out) :: data_out
622  sll_int32, dimension(:), intent(in) :: dorder
623 
627  sll_real64, dimension(:), intent(in) ::displacement
628 
629  sll_int32 :: i1, i2, k2, counter, j, index_parent, index_parent_old
630  sll_int32, dimension(4) :: no, ind_order
631  sll_int32, dimension(:, :), allocatable :: ind
632  sll_real64 :: coordinate_self, coordinate_ancestor, factor, disp
633 
634  sll_allocate(ind(interpolator%max_level + 1, 2), i1);
635  ! Dehierarchization along dimension dorder(1) only
636  do j = interpolator%order, 2, -1
637  call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part_order &
638  (data_in, 1, 1, dorder, j)
639  end do
640 
641  call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part(data_in, 1, 1, dorder)
642 
643  ! Interpolation in dorder(1)/dorder(2)-plane
644  ind_order(dorder(1)) = 0
645  no(dorder(1)) = 1
646  do i2 = 0, interpolator%levels(dorder(2))
647  ind_order(dorder(2)) = i2
648  no(dorder(2)) = max(2**(i2 - 1), 1);
649  ind(1, dorder(1)) = 0;
650  do k2 = 0, no(dorder(2)) - 1
651  ind(ind_order(dorder(2)) + 1, dorder(2)) = k2;
652  disp = displacement(2**(i2 - 1) + k2 + 1)
653 
654  counter = interpolator%index( &
655  ind_order(1), ind_order(2)) + &
656  ind(ind_order(1) + 1, 1)*no(2) &
657  + ind(ind_order(2) + 1, 2)
658 
659  ! Evaluate along dorder(1)-stripe
660  call interpolator%sll_t_sparse_grid_interpolator%interpolate_disp_1d_periodic_self &
661  (disp, dorder(1), &
662  min(interpolator%levels(dorder(1)), interpolator%max_level &
663  - ind_order(dorder(2))), counter, data_in, data_out)
664  ! Evaluate hierarchically along dorder(2) dimension (dorder(1)-stripe-wise)
665  index_parent = max( &
666  interpolator%hierarchy(counter)%parent(2*dorder(2) - 1), &
667  interpolator%hierarchy(counter)%parent(2*dorder(2)))
668  coordinate_self = interpolator%hierarchy(counter)%coordinate(dorder(2))
669  index_parent_old = counter;
670  do while (index_parent < index_parent_old)
671  coordinate_ancestor = interpolator%hierarchy(index_parent)% &
672  coordinate(dorder(2))
673  call interpolator%basis_function((coordinate_self - coordinate_ancestor)/ &
674  interpolator%length(dorder(2))* &
675  real(2**(interpolator%hierarchy(index_parent)%level(dorder(2))), f64), &
676  factor, &
677  interpolator%hierarchy(index_parent)%function_type(dorder(2)))
678  call interpolator%sll_t_sparse_grid_interpolator%interpolate_disp_1d_periodic_for_neighbor &
679  (disp, factor, &
680  dorder(1), min(interpolator%levels(dorder(1)), &
681  interpolator%max_level - &
682  interpolator%hierarchy(index_parent)%level(dorder(2))), &
683  min(interpolator%levels(dorder(1)), interpolator%max_level &
684  - ind_order(dorder(2))), index_parent, counter, &
685  data_in, data_out)
686  index_parent_old = index_parent;
687  index_parent = &
688  max(interpolator%hierarchy(index_parent)%parent(2*dorder(2) - 1), &
689  interpolator%hierarchy(index_parent)%parent(2*dorder(2)))
690  end do
691 
692  end do
693  end do
694 
695  end subroutine interpolate_disp_nconst_in_1d
696 
698  subroutine interpolate_disp_linnconst_in_1d(interpolator, displacement, dorder, data_in, data_out)
699  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
700  sll_real64, dimension(:), intent(inout) :: data_in
701  sll_real64, dimension(:), intent(out) :: data_out
702  sll_int32, dimension(:), intent(in) :: dorder
703  sll_real64, intent(in) ::displacement
704  sll_int32 :: i1, i2, k2, counter, j, index_parent, index_parent_old
705  sll_int32, dimension(4) :: no, ind_order
706  sll_int32, dimension(:, :), allocatable :: ind
707  sll_real64 :: coordinate_self, coordinate_ancestor, factor, disp
708 
709  sll_allocate(ind(interpolator%max_level + 1, 2), i1);
710  ! Dehierarchization along dimension dorder(1) only
711  do j = interpolator%order, 2, -1
712  call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part_order &
713  (data_in, 1, 1, dorder, j)
714  end do
715 
716  call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part(data_in, 1, 1, dorder)
717 
718  ! Interpolation in dorder(1)/dorder(2)-plane
719  ind_order(dorder(1)) = 0
720  no(dorder(1)) = 1
721  do i2 = 0, interpolator%levels(dorder(2))
722  ind_order(dorder(2)) = i2
723  no(dorder(2)) = max(2**(i2 - 1), 1);
724  ind(1, dorder(1)) = 0;
725  do k2 = 0, no(dorder(2)) - 1
726  ind(ind_order(dorder(2)) + 1, dorder(2)) = k2;
727  counter = interpolator%index( &
728  ind_order(1), ind_order(2)) + &
729  ind(ind_order(1) + 1, 1)*no(2) &
730  + ind(ind_order(2) + 1, 2)
731  disp = displacement*interpolator%hierarchy(counter)%coordinate(dorder(2));
732  ! Evaluate along dorder(1)-stripe
733  call interpolator%sll_t_sparse_grid_interpolator%interpolate_disp_1d_periodic_self &
734  (disp, dorder(1), &
735  min(interpolator%levels(dorder(1)), interpolator%max_level - &
736  ind_order(dorder(2))), counter, data_in, data_out)
737  ! Evaluate hierarchically along dorder(2) dimension (dorder(1)-stripe-wise)
738  index_parent = max( &
739  interpolator%hierarchy(counter)%parent(2*dorder(2) - 1), &
740  interpolator%hierarchy(counter)%parent(2*dorder(2)))
741  coordinate_self = interpolator%hierarchy(counter)%coordinate(dorder(2))
742  index_parent_old = counter;
743  do while (index_parent < index_parent_old)
744  coordinate_ancestor = interpolator%hierarchy(index_parent)% &
745  coordinate(dorder(2))
746  call interpolator%basis_function((coordinate_self - coordinate_ancestor)/ &
747  interpolator%length(dorder(2))* &
748  real(2**(interpolator%hierarchy(index_parent)%level(dorder(2))), f64), &
749  factor, &
750  interpolator%hierarchy(index_parent)%function_type(dorder(2)))
751  call interpolator%sll_t_sparse_grid_interpolator%interpolate_disp_1d_periodic_for_neighbor &
752  (disp, factor, &
753  dorder(1), min(interpolator%levels(dorder(1)), &
754  interpolator%max_level - &
755  interpolator%hierarchy(index_parent)%level(dorder(2))), &
756  min(interpolator%levels(dorder(2)), interpolator%max_level &
757  - ind_order(dorder(2))), index_parent, counter, &
758  data_in, data_out)
759  index_parent_old = index_parent;
760  index_parent = &
761  max(interpolator%hierarchy(index_parent)%parent(2*dorder(2) - 1), &
762  interpolator%hierarchy(index_parent)%parent(2*dorder(2)))
763  end do
764 
765  end do
766  end do
767 
768  end subroutine interpolate_disp_linnconst_in_1d
769 
770 !------------------------------------------------------------------------------!
771 ! Functions to evaluate fg on sg and sg on fg
772 
774  subroutine fg_to_sg(interpolator, fg_values, sg_values)
775  sll_real64, dimension(:, :), intent(in) :: fg_values
776  sll_real64, dimension(:), intent(out) :: sg_values
777  class(sll_t_sparse_grid_2d), intent(in) :: interpolator
778 
779  sll_int32 :: j
780  sll_int32, dimension(2) :: fg_ind
781 
782  do j = 1, interpolator%size_basis
783  fg_ind = fg_index(interpolator, j);
784  sg_values(j) = fg_values(fg_ind(1), fg_ind(2));
785  end do
786 
787  end subroutine fg_to_sg
788 
789 !PN DEFINED BUT NOT USED
790 !! Complex version of \a fg_to_sg
791 !subroutine sg_to_fg_complex(interpolator,sg_values,fg_values)
792 !sll_comp64, dimension(:,:), intent(out) :: fg_values
793 !sll_comp64, dimension(:), intent(in) :: sg_values
794 !class(sll_t_sparse_grid_2d), intent(in) :: interpolator
795 !sll_int32 :: j
796 !sll_int32, dimension(2) :: fg_ind
797 !
798 !do j=1,interpolator%size_basis
799 ! fg_ind = fg_index(interpolator,j);
800 ! fg_values(fg_ind(1),fg_ind(2)) = sg_values(j);
801 !end do
802 !
803 !end subroutine sg_to_fg_complex
804 
806  function fg_index(interpolator, sg_index)
807  sll_int32, intent(in) :: sg_index
808  sll_int32, dimension(2) :: fg_index
809  class(sll_t_sparse_grid_2d), intent(in) :: interpolator
810  sll_int32 :: j
811 
812  do j = 1, interpolator%dim
813  if (interpolator%hierarchy(sg_index)%level(j) == 0) then
814  fg_index(j) = 1;
815  else
816  fg_index(j) = 2**(interpolator%levels(j) - &
817  interpolator%hierarchy(sg_index)%level(j))* &
818  (1 + 2*interpolator%hierarchy(sg_index)%index_on_level(j)) + 1;
819  end if
820  end do
821 
822  end function fg_index
823 
824 ! End functions fg_to_sg and sg_to_fg
825 !------------------------------------------------------------------------------!
826 
827 !------------------------------------------------------------------------------!
828 !!!! SGFFT helper functions !!!!
829 
831  subroutine tohierarchical(interpolator, data_in, data_out)
832  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
833  sll_real64, dimension(:), intent(in) :: data_in
834  sll_comp64, dimension(:), intent(out) :: data_out
835  sll_int32 :: i, j
836 
837 ! Maybe possible with index in index to make dimension independent
838  do i = 0, interpolator%levels(2)
839  do j = interpolator%index(0, i), &
840  interpolator%index(0, i) + max(2**(i - 1), 1) - 1
841  call interpolator%ToHierarchical1D(1, &
842  min(interpolator%levels(1), interpolator%max_level - i), &
843  j, data_in, data_out)
844  end do
845  end do
846  do i = 0, interpolator%levels(1)
847  do j = interpolator%index(i, 0), &
848  interpolator%index(i, 0) + max(2**(i - 1), 1) - 1
849  call interpolator%ToHierarchical1D_comp(2, &
850  min(interpolator%levels(2), interpolator%max_level - i), j, data_out)
851  end do
852  end do
853 
854  end subroutine tohierarchical
855 
857  subroutine todehi(interpolator, data_array)
858  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
859  sll_comp64, dimension(:), intent(inout) :: data_array
860 
861  sll_int32 :: i, j
862 
863 ! Maybe possible with index in index to make dimension independent
864  do i = 0, interpolator%levels(1)
865  do j = interpolator%index(i, 0), &
866  interpolator%index(i, 0) + max(2**(i - 1), 1) - 1
867  call interpolator%ToDehi1D(2, &
868  min(interpolator%levels(2), interpolator%max_level - i), j, &
869  data_array)
870  end do
871  end do
872  do i = 0, interpolator%levels(2)
873  do j = interpolator%index(0, i), &
874  interpolator%index(0, i) + max(2**(i - 1), 1) - 1
875  call interpolator%ToDehi1D(1, &
876  min(interpolator%levels(1), interpolator%max_level - i), j, &
877  data_array)
878  end do
879  end do
880 
881  end subroutine todehi
882 
883  subroutine tohira(interpolator, data_array)
884  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
885  sll_comp64, dimension(:), intent(inout) :: data_array
886  sll_int32 :: i, j
887 
888  do i = 0, interpolator%levels(2)
889  do j = interpolator%index(0, i), &
890  interpolator%index(0, i) + max(2**(i - 1), 1) - 1
891  call interpolator%ToHira1D(1, &
892  min(interpolator%levels(1), interpolator%max_level - i), j, &
893  data_array)
894  end do
895  end do
896  do i = 0, interpolator%levels(1)
897  do j = interpolator%index(i, 0), &
898  interpolator%index(i, 0) + max(2**(i - 1), 1) - 1
899  call interpolator%ToHira1D(2, &
900  min(interpolator%levels(2), interpolator%max_level - i), j, &
901  data_array)
902  end do
903  end do
904 
905  end subroutine tohira
906 
907  subroutine tonodal(interpolator, data_in, data_out)
908  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
909  sll_comp64, dimension(:), intent(inout) :: data_in
910  sll_real64, dimension(:), intent(out) :: data_out
911  sll_int32 :: i, j
912 
913  do i = 0, interpolator%levels(1)
914  do j = interpolator%index(i, 0), &
915  interpolator%index(i, 0) + max(2**(i - 1), 1) - 1
916  call interpolator%ToNodal1D_comp(2, &
917  min(interpolator%levels(2), interpolator%max_level - i), &
918  j, data_in)
919  !call ToNodal1D(interpolator,2,interpolator%sparsegrid%levels-i,&
920  ! j,data_in,data_out)
921  end do
922  end do
923 
924  do i = 0, interpolator%levels(2)
925  do j = interpolator%index(0, i), &
926  interpolator%index(0, i) + max(2**(i - 1), 1) - 1
927  call interpolator%ToNodal1D(1, &
928  min(interpolator%levels(1), interpolator%max_level - i), &
929  j, data_in, data_out)
930  end do
931  end do
932 
933  end subroutine tonodal
934 
936  subroutine displace(interpolator, dim, displacement, data)
937  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
938  sll_comp64, dimension(:), intent(inout) :: data
939  sll_real64, intent(in) :: displacement
940  sll_int32, intent(in) :: dim
941 
942  sll_int32 :: i, j
943 
944  if (dim == 1) then
945  ! Maybe possible with index in index to make dimension independent
946  do i = 0, interpolator%levels(2)
947  do j = interpolator%index(0, i), &
948  interpolator%index(0, i) + max(2**(i - 1), 1) - 1
949  call interpolator%Displace1D(1, interpolator%levels(1) - i, &
950  j, displacement, data)
951  end do
952  end do
953  else
954  do i = 0, interpolator%levels(1)
955  do j = interpolator%index(i, 0), &
956  interpolator%index(i, 0) + max(2**(i - 1), 1) - 1
957  call interpolator%Displace1D(2, interpolator%levels(2) - i, &
958  j, displacement, data)
959  end do
960  end do
961  end if
962 
963  end subroutine displace
964 
965 !PN DEFINED BUT NOT USED
966 !subroutine DisplaceVar(interpolator,alpha1,alpha2,data)
967 ! class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
968 ! sll_comp64, dimension(:), intent(inout) :: data
969 ! sll_real64,dimension(:) ,intent(in) :: alpha1,alpha2
970 ! sll_int32 :: i,j,counter
971 !
972 !! Maybe possible with index in index to make dimension independent
973 ! counter = 1
974 ! do i = 0,interpolator%levels(2)
975 ! do j = interpolator%index(0,i),&
976 ! interpolator%index(0,i) + max(2**(i-1),1)-1
977 ! call interpolator%Displace1D(1,interpolator%levels(1)-i,&
978 ! j,alpha1(counter),data)
979 ! counter = counter +1
980 ! end do
981 ! end do
982 ! counter = 1
983 ! do i = 0,interpolator%levels(1)
984 ! do j = interpolator%index(i,0),&
985 ! interpolator%index(i,0) + max(2**(i-1),1)-1
986 ! call interpolator%Displace1D(2,interpolator%levels(2)-i,&
987 ! j,alpha2(counter),data)
988 ! counter = counter +1
989 ! end do
990 ! end do
991 !
992 !
993 !end subroutine DISPLACEVAR
994 
996  subroutine spfft(interpolator, data_in, data_out)
997  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
998  sll_real64, dimension(:), intent(in) :: data_in
999  sll_comp64, dimension(:), intent(out) :: data_out
1000 
1001  call interpolator%ToHierarchical(data_in, data_out)
1002  call interpolator%ToDehi(data_out)
1003 
1004  end subroutine spfft
1005 
1007  subroutine ispfft(interpolator, data_in, data_out)
1008  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
1009  sll_comp64, dimension(:), intent(inout) :: data_in
1010  sll_real64, dimension(:), intent(out) :: data_out
1011 
1012  call interpolator%toHira(data_in)
1013  call interpolator%ToNodal(data_in, data_out)
1014 
1015  end subroutine ispfft
1016 
1017 !!!! End SGFFT helper functions !!!!
1018 !------------------------------------------------------------------------------!
1019 
1020 !------------------------------------------------------------------------------!
1021 !!!! Filter functions !!!!
1022 
1023  subroutine filter_highest(interpolator, data)
1024  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
1025  sll_real64, dimension(:), intent(inout) :: data
1026  sll_int32 :: lev, l, j, ind
1027 
1028  lev = interpolator%levels(1);
1029  do l = 0, interpolator%max_level - lev
1030  ind = interpolator%index(0, l)
1031  call interpolator%extract_periodic(1, lev, ind, data, interpolator%stripe)
1032  call interpolator%hierarchical_stripe(interpolator%stripe, lev);
1033  do j = 2, 2**lev, 2
1034  interpolator%stripe(j) = interpolator%stripe(j)/2;
1035  end do
1036  call interpolator%dehierarchical_stripe(interpolator%stripe, lev);
1037  call interpolator%insert_periodic(1, lev, ind, interpolator%stripe, data);
1038  end do
1039 
1040  lev = interpolator%levels(2);
1041  do l = 0, interpolator%max_level - lev
1042  ind = interpolator%index(l, 0)
1043  call interpolator%extract_periodic(2, lev, ind, data, interpolator%stripe)
1044  call interpolator%hierarchical_stripe(interpolator%stripe, lev);
1045  do j = 2, 2**lev, 2
1046  interpolator%stripe(j) = interpolator%stripe(j)*0.5_f64;
1047  end do
1048  call interpolator%dehierarchical_stripe(interpolator%stripe, lev);
1049  call interpolator%insert_periodic(2, lev, ind, interpolator%stripe, data);
1050  end do
1051 
1052  end subroutine filter_highest
1053 
1054  subroutine filter(interpolator, data)
1055  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
1056  sll_real64, dimension(:), intent(inout) :: data
1057  sll_int32 :: j
1058 
1059  call interpolator%compute_hierarchical_surplus(data);
1060  do j = interpolator%level_mapping(interpolator%max_level) + 1, interpolator%size_basis, 2
1061  data(j) = data(j)*0.5_f64;
1062  end do
1063  call interpolator%compute_dehierarchical(data);
1064  end subroutine filter
1065 
1066 !PN DEFINED BUT NOT USED
1067 !subroutine filter_oscillations(interpolator, data)
1068 ! class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
1069 ! sll_real64, dimension(:), intent(inout) :: data
1070 ! sll_int32 :: j
1071 !
1072 ! call interpolator%compute_hierarchical_surplus(data);
1073 ! do j=interpolator%level_mapping(interpolator%max_level)+1, interpolator%size_basis,2
1074 ! data(j) = data(j)*0.5_f64;
1075 ! end do
1076 ! call interpolator%compute_dehierarchical(data);
1077 !
1078 !end subroutine filter_oscillations
1079 
1080  subroutine linear_filter(interpolator, data, hs, width)
1081  class(sll_t_sparse_grid_2d), intent(inout) :: interpolator
1082  sll_real64, dimension(:), intent(inout) :: data, hs
1083  sll_real64, dimension(2), intent(in) :: width
1084  sll_real64, dimension(2) :: dx, dxn
1085  sll_int32 :: j
1086 
1087  do j = 1, interpolator%size_basis
1088  dx = interpolator%hierarchy(j)%coordinate
1089  dxn(1) = dx(1); dxn(2) = dx(2) + width(2)
1090  data(j) = data(j)*0.5_f64 + 0.125_f64*interpolator%interpolate_from_interpolant_value(hs, dxn);
1091  dxn(2) = dx(2) - width(2)
1092  data(j) = data(j) + 0.125_f64*interpolator%interpolate_from_interpolant_value(hs, dxn);
1093  dxn(1) = dx(1) + width(1); dxn(2) = dx(2)
1094  data(j) = data(j) + 0.125_f64*interpolator%interpolate_from_interpolant_value(hs, dxn);
1095  dxn(1) = dx(1) - width(1)
1096  data(j) = data(j) + 0.125_f64*interpolator%interpolate_from_interpolant_value(hs, dxn);
1097  end do
1098 
1099  end subroutine linear_filter
1100 
1101 !!!! End filter functions !!!!
1102 !------------------------------------------------------------------------------!
1103 
1104 end module sll_m_sparse_grid_2d
Fortran module where set some physical and mathematical constants.
real(kind=f64), parameter, public sll_p_pi
Implementation of a 2D sparse grid with interpolation routines.
subroutine displace(interpolator, dim, displacement, data)
Compute the Fourier coefficients of at displaced grid points from Fourier coefficients.
subroutine interpolate_disp_nconst_in_1d(interpolator, displacement, dorder, data_in, data_out)
Functionality: Interpolates the function values for a displacement on in dimension (periodic b....
subroutine set_hierarchy_info(interpolator, counter, cdim, lvecin, kvecin, novecin)
Helfer function for initialization. Setting all the information needed for node counter of the sparse...
subroutine spfft(interpolator, data_in, data_out)
Fourier transform on sparse grid.
subroutine interpolate_array_disp_sgfft(interpolator, dim, displacment_in, data_in, data_out)
Compute value at displaced grid points using trigonometric interpolation (based on SG FFT)
subroutine fg_to_sg(interpolator, fg_values, sg_values)
Set sparse grid values from fg vector.
subroutine ispfft(interpolator, data_in, data_out)
Inverse Fourier transform.
subroutine todehi(interpolator, data_array)
subroutine interpolate_const_disp(interpolator, dorder, displacement, data_in, data_out, hiera)
Interpolation function for interpolation at (constantly) displaced grid points; displacement only in ...
subroutine filter_highest(interpolator, data)
subroutine filter(interpolator, data)
subroutine tohierarchical(interpolator, data_in, data_out)
Compute Fourier coefficient on sparse grid.
subroutine tohira(interpolator, data_array)
integer(kind=i32) function, dimension(2) fg_index(interpolator, sg_index)
Compute the index of a sparse grid node on level "level" with index "index_on_level" on full grid wit...
real(kind=f64) function interpolate_value_sg(interpolator, data, eta)
Value at eta interpolated from the hierarchical surplus data using standard sparse grid interpolation...
subroutine initialize_sg2d(interpolator, levels, order, interpolation, interpolation_type, eta_min, eta_max, boundary, modified)
Initialize a 2d sparse grid object.
subroutine set_hierarchy_info_boundary(interpolator, counter, cdim, lvecin, kvecin, novecin)
Same as set_hierarchy_info but for points on the boundary along dimension cdim.
real(kind=f64) function interpolate_from_hierarchical_surplus(interpolator, data, eta)
Implementation of interpolate_value_sg for periodic sparse grid.
subroutine tonodal(interpolator, data_in, data_out)
subroutine linear_filter(interpolator, data, hs, width)
subroutine interpolate_disp_linnconst_in_1d(interpolator, displacement, dorder, data_in, data_out)
As interpolate_disp_nconst_in_1d but displacement dependent on displacement*coordinate(dorder(2))
real(kind=f64) function interpolate_from_hierarchical_surplus_boundary(interpolator, data, eta)
Implementation of interpolate_value_sg for sparse grid with boundary.
Dimension-independent functions for sparse grid with polynomial basis functions.
Sparse grid object for 2d with interpolation routines.
    Report Typos and Errors