Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_sparse_grid_4d.F90
Go to the documentation of this file.
1 
5 
6 module sll_m_sparse_grid_4d
7 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8 #include "sll_memory.h"
9 #include "sll_working_precision.h"
10 
11  use sll_m_constants, only: &
12  sll_p_pi
13 
14  use sll_m_sparse_grid_interpolator, only: &
15  sll_t_sparse_grid_interpolator
16 
17  implicit none
18 
19  public :: &
20  sll_t_sparse_grid_4d
21 
22  private
23 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
24 
26  type, extends(sll_t_sparse_grid_interpolator) :: sll_t_sparse_grid_4d
27  sll_int32, dimension(:, :, :, :), pointer :: index
28 
29  contains
30  procedure :: init => initialize_sg4d! Initialization routine
31  procedure :: interpolate_from_interpolant_value ! Compute the value of the sparse grid interpolant at position eta
32  procedure :: interpolate_disp_nconst_in_1d ! Interpolate along one (x)-direction with displacement non-constant in one (v)-direction
33  procedure :: interpolate_disp_linnconst_in_1d => interpolate4d_disp_linnconst_in_1d
34  procedure :: interpolate_disp_nconst_in_2d ! Interpolate along one (v)-direction with displacement non-constant in all x-directions
35  procedure :: interpolate_const_disp
36 
37  end type sll_t_sparse_grid_4d
38 
39 contains
40 
41 !------------------------------------------------------------------------------!
42 !!!! Interpolation routines !!!!
43 
45  function interpolate_from_interpolant_value(interpolator, data, eta) result(val)
46  class(sll_t_sparse_grid_4d), intent(inout) :: interpolator
47  sll_real64 :: val
48  sll_real64, dimension(:), intent(in) :: data
49  sll_real64, dimension(:), intent(in) :: eta
50  val = interpolate_from_hierarchical_surplus( &
51  interpolator, data, eta)
52 
53  end function interpolate_from_interpolant_value
54 
60  subroutine interpolate_const_disp(interpolator, dorder, displacement, data_in, data_out, hiera)
61  class(sll_t_sparse_grid_4d), intent(inout) :: interpolator
62  sll_real64, dimension(:), intent(inout) :: data_in
63  sll_real64, dimension(:), intent(out) :: data_out
64  sll_int32, dimension(:), intent(in) :: dorder
65  sll_real64, intent(in) ::displacement
66  logical, intent(in) :: hiera
67 
68  sll_int32 :: i1, i2, i3, i4, k2, k3, k4, counter, j
69  sll_int32, dimension(4) :: l, no, ind_order
70  sll_int32, dimension(:, :), allocatable :: ind
71 
72  sll_allocate(ind(interpolator%max_level + 1, 4), i1);
73  ind_order(dorder(1)) = 0
74  no(dorder(1)) = 1
75  do i3 = 0, interpolator%levels(dorder(3))
76  ind_order(dorder(3)) = i3
77  l(dorder(3)) = i3
78  no(dorder(3)) = max(2**(i3 - 1), 1);
79  do i4 = 0, min(interpolator%max_level - i3, interpolator%levels(dorder(4)))
80  ind_order(dorder(4)) = i4
81  l(dorder(4)) = i4
82  no(dorder(4)) = max(2**(i4 - 1), 1);
83  do i2 = 0, min(interpolator%max_level - i3 - i4, interpolator%levels(dorder(2)))
84  ind_order(dorder(2)) = i2
85  no(dorder(2)) = max(2**(i2 - 1), 1);
86  ind(1, dorder(1)) = 0;
87  do k2 = 0, no(dorder(2)) - 1
88  ind(ind_order(dorder(2)) + 1, dorder(2)) = k2;
89  do k3 = 0, no(dorder(3)) - 1
90  ind(ind_order(dorder(3)) + 1, dorder(3)) = k3;
91  do k4 = 0, no(dorder(4)) - 1
92  ind(ind_order(dorder(4)) + 1, dorder(4)) = k4;
93  counter = interpolator%index( &
94  ind_order(1), ind_order(2), ind_order(3), ind_order(4)) + &
95  ind(ind_order(1) + 1, 1)*no(2)*no(3)*no(4) &
96  + ind(ind_order(2) + 1, 2)*no(3)*no(4) + &
97  ind(ind_order(3) + 1, 3)*no(4) + ind(ind_order(4) + 1, 4)
98  ! Evaluate along dorder(1)-stripe
99  call interpolator%interpolate_disp_1d_periodic &
100  (displacement, dorder(1), &
101  min(interpolator%levels(dorder(1)), &
102  interpolator%max_level - ind_order(dorder(2)) - &
103  ind_order(dorder(3)) - &
104  ind_order(dorder(4))), counter, data_in, data_out, hiera)
105  end do
106  end do
107  end do
108  end do
109  end do
110  end do
111 
112  if (hiera .EQV. .false.) then
113  ! Dehierarchization along dimension dorder(1) only
114  do j = interpolator%order, 2, -1
115  call interpolator%dehierarchical_part_order &
116  (data_out, &
117  interpolator%dim, 2, dorder, j)
118  end do
119 
120  call interpolator%dehierarchical_part(data_out, &
121  interpolator%dim, 2, dorder)
122  end if
123 
124  end subroutine interpolate_const_disp
125 
130  subroutine interpolate_disp_nconst_in_1d(interpolator, displacement, dorder, data_in, data_out)
131  class(sll_t_sparse_grid_4d), intent(inout) :: interpolator
132  sll_real64, dimension(:), intent(inout) :: data_in
133  sll_real64, dimension(:), intent(out) :: data_out
134  sll_int32, dimension(:), intent(in) :: dorder
135 
136  sll_real64, dimension(:), intent(in) ::displacement
137  sll_int32 :: i1, i2, i3, i4, k2, k3, k4, counter, j, index_parent, index_parent_old
138  sll_int32, dimension(4) :: l, no, ind_order
139  sll_int32, dimension(:, :), allocatable :: ind
140  sll_real64 :: coordinate_self, coordinate_ancestor, factor, disp
141 
142  sll_allocate(ind(interpolator%max_level + 1, 4), i1);
143  ! Dehierarchization along dimension dorder(1) only
144  do j = interpolator%order, 2, -1
145  call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part_order &
146  (data_in, 1, 1, dorder, j)
147  end do
148 
149  call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part(data_in, 1, 1, dorder)
150  ! Interpolation in dorder(1)/dorder(2)-plane
151  ind_order(dorder(1)) = 0
152  no(dorder(1)) = 1
153  do i3 = 0, interpolator%levels(dorder(3))
154  ind_order(dorder(3)) = i3
155  l(dorder(3)) = i3
156  no(dorder(3)) = max(2**(i3 - 1), 1);
157  do i4 = 0, min(interpolator%max_level - i3, interpolator%levels(dorder(4)))
158  ind_order(dorder(4)) = i4
159  l(dorder(4)) = i4
160  no(dorder(4)) = max(2**(i4 - 1), 1);
161  do i2 = 0, min(interpolator%max_level - i3 - i4, interpolator%levels(dorder(2)))
162  ind_order(dorder(2)) = i2
163  no(dorder(2)) = max(2**(i2 - 1), 1);
164  ind(1, dorder(1)) = 0;
165  do k2 = 0, no(dorder(2)) - 1
166  ind(ind_order(dorder(2)) + 1, dorder(2)) = k2;
167  disp = displacement(2**(i2 - 1) + k2 + 1)
168  do k3 = 0, no(dorder(3)) - 1
169  ind(ind_order(dorder(3)) + 1, dorder(3)) = k3;
170  do k4 = 0, no(dorder(4)) - 1
171  ind(ind_order(dorder(4)) + 1, dorder(4)) = k4;
172  counter = interpolator%index( &
173  ind_order(1), ind_order(2), ind_order(3), ind_order(4)) + &
174  ind(ind_order(1) + 1, 1)*no(2)*no(3)*no(4) &
175  + ind(ind_order(2) + 1, 2)*no(3)*no(4) + &
176  ind(ind_order(3) + 1, 3)*no(4) + ind(ind_order(4) + 1, 4)
177  ! Evaluate along dorder(1)-stripe
178  call interpolator%sll_t_sparse_grid_interpolator%interpolate_disp_1d_periodic_self &
179  (disp, dorder(1), &
180  min(interpolator%levels(dorder(1)), &
181  interpolator%max_level - ind_order(dorder(2)) - &
182  ind_order(dorder(3)) - &
183  ind_order(dorder(4))), counter, data_in, data_out)
184  ! Evaluate hierarchically along dorder(2) dimension (dorder(1)-stripe-wise)
185  index_parent = max(interpolator%hierarchy(counter)%parent(2*dorder(2) - 1), &
186  interpolator%hierarchy(counter)%parent(2*dorder(2)))
187  coordinate_self = interpolator%hierarchy(counter)%coordinate(dorder(2))
188  index_parent_old = counter;
189  do while (index_parent < index_parent_old)
190  coordinate_ancestor = interpolator%hierarchy(index_parent)% &
191  coordinate(dorder(2))
192  call interpolator%basis_function((coordinate_self - coordinate_ancestor)/ &
193  interpolator%length(dorder(2))* &
194  real(2**(interpolator%hierarchy(index_parent)%level(dorder(2))), f64), factor, &
195  interpolator%hierarchy(index_parent)%function_type(dorder(2)))
196  call interpolator%sll_t_sparse_grid_interpolator%interpolate_disp_1d_periodic_for_neighbor &
197  (disp, factor, &
198  dorder(1), min(interpolator%levels(dorder(1)), &
199  interpolator%max_level - &
200  interpolator%hierarchy(index_parent)%level(dorder(2)) - &
201  ind_order(dorder(3)) - ind_order(dorder(4))), &
202  min(interpolator%levels(dorder(1)), &
203  interpolator%max_level - ind_order(dorder(2)) - &
204  ind_order(dorder(3)) - &
205  ind_order(dorder(4))), index_parent, counter, &
206  data_in, data_out)
207  index_parent_old = index_parent;
208  index_parent = &
209  max(interpolator%hierarchy(index_parent)%parent(2*dorder(2) - 1), &
210  interpolator%hierarchy(index_parent)%parent(2*dorder(2)))
211  end do
212 
213  end do
214  end do
215  end do
216  end do
217  end do
218  end do
219 
220  ! Dehierarchization along dimension dorder(3) dorder(4)
221  do j = interpolator%order, 2, -1
222  call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part_order &
223  (data_out, 4, 3, dorder, j)
224  end do
225 
226  call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part(data_out, 4, 3, dorder)
227 
228  end subroutine interpolate_disp_nconst_in_1d
229 
231  subroutine interpolate4d_disp_linnconst_in_1d(interpolator, displacement, dorder, data_in, data_out)
232  class(sll_t_sparse_grid_4d), intent(inout) :: interpolator
233  sll_real64, dimension(:), intent(inout) :: data_in
234  sll_real64, dimension(:), intent(out) :: data_out
235  sll_int32, dimension(:), intent(in) :: dorder
236  sll_real64, intent(in) ::displacement
237  sll_int32 :: i1, i2, i3, i4, k2, k3, k4, counter, j, index_parent, index_parent_old
238  sll_int32, dimension(4) :: l, no, ind_order
239  sll_int32, dimension(:, :), allocatable :: ind
240  sll_real64 :: coordinate_self, coordinate_ancestor, factor, disp
241 
242  sll_allocate(ind(interpolator%max_level + 1, 4), i1);
243  ! Dehierarchization along dimension dorder(1) only
244  do j = interpolator%order, 2, -1
245  call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part_order &
246  (data_in, 1, 1, dorder, j)
247  end do
248 
249  call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part(data_in, 1, 1, dorder)
250 
251  ! Interpolation in dorder(1)/dorder(2)-plane
252  ind_order(dorder(1)) = 0
253  no(dorder(1)) = 1
254  do i3 = 0, interpolator%levels(dorder(3))
255  ind_order(dorder(3)) = i3
256  l(dorder(3)) = i3
257  no(dorder(3)) = max(2**(i3 - 1), 1);
258  do i4 = 0, min(interpolator%max_level - i3, interpolator%levels(dorder(4)))
259  ind_order(dorder(4)) = i4
260  l(dorder(4)) = i4
261  no(dorder(4)) = max(2**(i4 - 1), 1);
262  do i2 = 0, min(interpolator%max_level - i3 - i4, interpolator%levels(dorder(2)))
263  ind_order(dorder(2)) = i2
264  no(dorder(2)) = max(2**(i2 - 1), 1);
265  ind(1, dorder(1)) = 0;
266  do k2 = 0, no(dorder(2)) - 1
267  ind(ind_order(dorder(2)) + 1, dorder(2)) = k2;
268  do k3 = 0, no(dorder(3)) - 1
269  ind(ind_order(dorder(3)) + 1, dorder(3)) = k3;
270  do k4 = 0, no(dorder(4)) - 1
271  ind(ind_order(dorder(4)) + 1, dorder(4)) = k4;
272  counter = interpolator%index( &
273  ind_order(1), ind_order(2), ind_order(3), ind_order(4)) + &
274  ind(ind_order(1) + 1, 1)*no(2)*no(3)*no(4) &
275  + ind(ind_order(2) + 1, 2)*no(3)*no(4) + &
276  ind(ind_order(3) + 1, 3)*no(4) + ind(ind_order(4) + 1, 4)
277  disp = displacement*interpolator%hierarchy(counter)%coordinate(dorder(2))
278  ! Evaluate along dorder(1)-stripe
279  call interpolator%sll_t_sparse_grid_interpolator%interpolate_disp_1d_periodic_self &
280  (disp, dorder(1), &
281  min(interpolator%levels(dorder(1)), interpolator%max_level &
282  - ind_order(dorder(2)) - ind_order(dorder(3)) - &
283  ind_order(dorder(4))), counter, data_in, data_out)
284  ! Evaluate hierarchically along dorder(2) dimension (dorder(1)-stripe-wise)
285  index_parent = max(interpolator%hierarchy(counter)%parent(2*dorder(2) - 1), &
286  interpolator%hierarchy(counter)%parent(2*dorder(2)))
287  coordinate_self = interpolator%hierarchy(counter)%coordinate(dorder(2))
288  index_parent_old = counter;
289  do while (index_parent < index_parent_old)
290  coordinate_ancestor = interpolator%hierarchy(index_parent)% &
291  coordinate(dorder(2))
292  call interpolator%basis_function((coordinate_self - coordinate_ancestor)/ &
293  interpolator%length(dorder(2))* &
294  real(2**(interpolator%hierarchy(index_parent)%level(dorder(2))), f64), factor, &
295  interpolator%hierarchy(index_parent)%function_type(dorder(2)))
296  call interpolator%sll_t_sparse_grid_interpolator%interpolate_disp_1d_periodic_for_neighbor &
297  (disp, factor, &
298  dorder(1), min(interpolator%levels(dorder(1)), &
299  interpolator%max_level - &
300  interpolator%hierarchy(index_parent)%level(dorder(2)) - &
301  ind_order(dorder(3)) - ind_order(dorder(4))), &
302  min(interpolator%levels(dorder(1)), &
303  interpolator%max_level - &
304  ind_order(dorder(2)) - ind_order(dorder(3)) - &
305  ind_order(dorder(4))), index_parent, counter, &
306  data_in, data_out)
307  index_parent_old = index_parent;
308  index_parent = &
309  max(interpolator%hierarchy(index_parent)%parent(2*dorder(2) - 1), &
310  interpolator%hierarchy(index_parent)%parent(2*dorder(2)))
311  end do
312 
313  end do
314  end do
315  end do
316  end do
317  end do
318  end do
319 
320  ! Dehierarchization along dimension dorder(3) dorder(4)
321  do j = interpolator%order, 2, -1
322  call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part_order &
323  (data_out, 4, 3, dorder, j)
324  end do
325 
326  call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part(data_out, 4, 3, dorder)
327 
328  end subroutine interpolate4d_disp_linnconst_in_1d
329 
331  subroutine interpolate_disp_nconst_in_2d(interpolator, displacement, dorder, data_in, data_out)
332  class(sll_t_sparse_grid_4d), intent(inout) :: interpolator
333  sll_real64, dimension(:), intent(inout) :: data_in
334  sll_real64, dimension(:), intent(out) :: data_out
335  sll_int32, dimension(:), intent(in) :: dorder
336  sll_real64, dimension(:), intent(in) ::displacement
337  sll_int32 :: i1, i2, i3, i4, k2, k3, k4, counter, j, index_parent, index_old, index_upper, index_upper_old
338  sll_int32 :: counter_disp, l23
339  sll_int32, dimension(4) :: l, no, ind_order
340  sll_int32, dimension(:, :), allocatable :: ind
341  sll_real64 :: coordinate_self, coordinate_ancestor, coordinate_self_upper, factor_upper, factor_lower, factor, disp
342 
343  sll_allocate(ind(interpolator%max_level + 1, 4), i1);
344  ! Dehierarchization along dimension dorder(1) only
345  do j = interpolator%order, 2, -1
346  call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part_order &
347  (data_in, 1, 1, dorder, j)
348  end do
349 
350  call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_part( &
351  data_in, 1, 1, dorder)
352 
353  !print*, data_in
354 
355  ! Interpolation in dorder(1)/dorder(2)-plane
356  ind_order(dorder(1)) = 0
357  no(dorder(1)) = 1
358  do i4 = 0, interpolator%levels(dorder(4))
359  ind_order(dorder(4)) = i4
360  l(dorder(4)) = i4
361  no(dorder(4)) = max(2**(i4 - 1), 1);
362  counter_disp = 0;
363  do l23 = 0, min(interpolator%max_level - i4, interpolator%levels(dorder(2)))
364  do i2 = 0, l23
365  ind_order(dorder(2)) = i2
366  no(dorder(2)) = max(2**(i2 - 1), 1);
367  i3 = l23 - i2
368  ind_order(dorder(3)) = i3
369  l(dorder(3)) = i3
370  no(dorder(3)) = max(2**(i3 - 1), 1);
371  ind(1, dorder(1)) = 0;
372  do k2 = 0, max(2**(ind_order(dorder(2)) - 1), 1) - 1
373  ind(ind_order(dorder(2)) + 1, dorder(2)) = k2;
374  do k3 = 0, max(2**(ind_order(dorder(3)) - 1), 1) - 1
375  counter_disp = counter_disp + 1;
376  disp = displacement(counter_disp);
377  ind(ind_order(dorder(3)) + 1, dorder(3)) = k3;
378  do k4 = 0, max(2**(ind_order(dorder(4)) - 1), 1) - 1
379  ind(ind_order(dorder(4)) + 1, dorder(4)) = k4;
380  ! print*, ind_order
381  ! print*, k2,k3,k4
382  ! print*, ind(ind_order(1)+1,1),ind(ind_order(2)+1,2),ind(ind_order(3)+1,3),ind(ind_order(4)+1,4)
383  ! print*, '------------------------------------------------- '
384  counter = interpolator%index( &
385  ind_order(1), ind_order(2), ind_order(3), ind_order(4)) + &
386  ind(ind_order(1) + 1, 1)*no(2)*no(3)*no(4) &
387  + ind(ind_order(2) + 1, 2)*no(3)*no(4) + &
388  ind(ind_order(3) + 1, 3)*no(4) + ind(ind_order(4) + 1, 4)
389 
390  ! Evaluate along dorder(1)-stripe
391  call interpolator%sll_t_sparse_grid_interpolator%interpolate_disp_1d_periodic_self &
392  (disp, dorder(1), &
393  min(interpolator%levels(dorder(1)), &
394  interpolator%max_level - ind_order(dorder(2)) - &
395  ind_order(dorder(3)) - &
396  ind_order(dorder(4))), counter, data_in, data_out)
397  ! Evaluate hierarchically along dorder(2)&dorder(3) dimension (dorder(1)-stripe-wise)
398  index_upper = counter
399  index_upper_old = index_upper + 1
400 
401  !print*, counter, counter_disp , disp
402  !print*, counter , interpolator%index(&
403  ! ind_order(1),ind_order(2),ind_order(3),ind_order(4)), ind_order
404 
405  coordinate_self = interpolator%hierarchy(counter)%coordinate(dorder(2))
406  coordinate_self_upper = interpolator%hierarchy(counter)%coordinate(dorder(3))
407  !print*, counter
408  do while (index_upper < index_upper_old)
409  coordinate_ancestor = interpolator%hierarchy(index_upper)% &
410  coordinate(dorder(3))
411  call interpolator%basis_function((coordinate_self_upper - coordinate_ancestor)/ &
412  interpolator%length(dorder(3))* &
413  2**(interpolator%hierarchy(index_upper)%level(dorder(3))), &
414  factor_upper, &
415  interpolator%hierarchy(index_upper)%function_type(dorder(3)))
416 
417  if (index_upper == counter) then
418  index_parent = &
419  max(interpolator%hierarchy(index_upper)%parent(2*dorder(2) - 1), &
420  interpolator%hierarchy(index_upper)%parent(2*dorder(2)))
421  index_old = index_upper
422  else
423  index_parent = index_upper;
424  index_old = index_parent + 1
425  end if
426 
427  do while (index_parent < index_old)
428  coordinate_ancestor = interpolator%hierarchy(index_parent)% &
429  coordinate(dorder(2))
430  call interpolator%basis_function((coordinate_self - coordinate_ancestor)/ &
431  interpolator%length(dorder(2))* &
432  real(2**(interpolator%hierarchy(index_parent)%level(dorder(2))), f64), &
433  factor_lower, &
434  interpolator%hierarchy(index_parent)%function_type(dorder(2)))
435  factor = factor_upper*factor_lower
436  !if(counter==30) then
437  ! print*, index_parent,data_out(30)
438  !end if
439  ! print*, counter, index_parent, factor, factor_upper, factor_lower
440  call interpolator%sll_t_sparse_grid_interpolator%interpolate_disp_1d_periodic_for_neighbor &
441  (disp, factor, &
442  dorder(1), min(interpolator%levels(dorder(1)), &
443  interpolator%max_level - &
444  interpolator%hierarchy(index_parent)%level(dorder(2)) - &
445  interpolator%hierarchy(index_parent)%level(dorder(3)) - &
446  ind_order(dorder(4))), &
447  min(interpolator%levels(dorder(1)), &
448  interpolator%max_level - ind_order(dorder(2)) - &
449  ind_order(dorder(3)) - &
450  ind_order(dorder(4))), index_parent, counter, &
451  data_in, data_out)
452  index_old = index_parent
453  index_parent = &
454  max(interpolator%hierarchy(index_parent)%parent(2*dorder(2) - 1), &
455  interpolator%hierarchy(index_parent)%parent(2*dorder(2)))
456  end do
457  index_upper_old = index_upper;
458  index_upper = &
459  max(interpolator%hierarchy(index_upper)%parent(2*dorder(3) - 1), &
460  interpolator%hierarchy(index_upper)%parent(2*dorder(3)));
461  end do
462  end do
463  end do
464  end do
465  end do
466  end do
467  end do
468 
469 ! Hierarchization along dimension dorder(1)-dorder(3)
470  call interpolator%sll_t_sparse_grid_interpolator%hierarchical_part(data_out, 3, 1, dorder)
471 
472  do j = 2, interpolator%order
473  call interpolator%sll_t_sparse_grid_interpolator%hierarchical_part_order &
474  (data_out, 3, 1, dorder, j)
475  end do
476 
477 ! Dehierarchization of the hierarchical surplus
478  do j = interpolator%order, 2, -1
479  call interpolator%sll_t_sparse_grid_interpolator%dehierarchical_order &
480  (data_out, j)
481  end do
482 
483  call interpolator%sll_t_sparse_grid_interpolator%dehierarchical &
484  (data_out)
485 
486  end subroutine interpolate_disp_nconst_in_2d
487 
488 !PN DEFINED BUT NOT USED
489 !! Note dorder should contain: dorder(1) dimension with displacement (1 or 2), dorder(2) dimension where displacement is non-constant (3 or 4), dorder(3) other of 1 or 2, dorder(4) other of 3 or 4
490 !subroutine displace_disp_nconst1(interpolator, surplus, data, dorder,displacement)
491 ! class(sll_t_sparse_grid_4d), intent(inout) :: interpolator
492 ! sll_real64, dimension(:), intent(in) :: surplus
493 ! sll_real64, dimension(:), intent(out) :: data
494 ! sll_int32, dimension(:), intent(in) :: dorder
495 ! sll_real64, dimension(:), intent(in) :: displacement
496 ! sll_int32 :: i3, i4, i2, i1, k1,k2,k3,k4, counter
497 ! sll_int32, dimension(4) :: ind_order,l
498 ! sll_real64, dimension(2) :: eta
499 ! sll_int32, dimension(4) :: no
500 ! sll_int32, dimension(:,:), allocatable :: ind
501 !
502 ! SLL_ALLOCATE(ind(interpolator%max_level,4), i1);
503 !
504 !
505 ! ind_order(dorder(1)) = 0
506 ! do i3 = 0,interpolator%levels(dorder(3))
507 ! ind_order(dorder(3)) = i3
508 ! l(dorder(3)) = i3
509 ! no(dorder(3)) = max(2**(i3-1),1);
510 ! do i4 = 0,min(interpolator%max_level - i3,interpolator%levels(dorder(4)))
511 ! ind_order(dorder(4)) = i4
512 ! l(dorder(4)) = i4
513 ! no(dorder(4)) = max(2**(i4-1),1);
514 ! do i2 = 0,min(interpolator%max_level - i3 -i4,interpolator%levels(dorder(2)))
515 ! ind_order(dorder(2)) = i2
516 ! no(dorder(2)) = max(2**(i2-1),1);
517 ! do i1 = 0,min(interpolator%max_level - i3 - i4 -i2,interpolator%levels(dorder(1)))
518 ! ind_order(dorder(1)) = i1
519 ! no(dorder(1)) = max(2**(i1-1),1);
520 ! do k1 = 0,max(2**(ind_order(dorder(1))-1),1)-1
521 ! ind(ind_order(dorder(1))+1,dorder(1)) = k1;
522 ! do k2 = 0,max(2**(ind_order(dorder(2))-1),1)-1
523 ! ind(ind_order(dorder(2))+1,dorder(2)) = k2;
524 ! do k3 = 0,max(2**(ind_order(dorder(3))-1),1)-1
525 ! ind(ind_order(dorder(3))+1,dorder(3)) = k3;
526 ! do k4 = 0,max(2**(ind_order(dorder(4))-1),1)-1
527 ! ind(ind_order(dorder(4))+1,dorder(4)) = k4;
528 ! counter = interpolator%index(&
529 ! ind_order(1),ind_order(2),ind_order(3),ind_order(4))+&
530 ! ind(ind_order(1)+1,1)*no(2)*no(3)*no(4)&
531 ! +ind(ind_order(2)+1,2)*no(3)*no(4)+&
532 ! ind(ind_order(3)+1,3)*no(4)+ind(ind_order(4)+1,4)
533 ! eta(1) = interpolator%hierarchy(counter)%coordinate(dorder(1))+&
534 ! displacement(2**(i2-1)+k2)
535 ! eta(2) = interpolator%hierarchy(counter)%coordinate(dorder(2))
536 !
537 ! l(dorder(2)) = ind_order(dorder(2))
538 ! data(counter) = interpolate_from_2D_hierarchical_surplus( &
539 ! interpolator,surplus, eta, dorder, no, ind, l );
540 ! end do
541 ! end do
542 ! end do
543 ! end do
544 ! end do
545 ! end do
546 ! end do
547 ! end do
548 !
549 ! call interpolator%sll_t_sparse_grid_interpolator%hierarchical_part(data,2,1,dorder)
550 !
551 !
552 ! call interpolator%sll_t_sparse_grid_interpolator%dehierarchical(data)
553 !
554 !
555 !end subroutine displace_disp_nconst1
556 
557 ! helper functions
558 
559  function interpolate_from_hierarchical_surplus(interpolator, data, eta) result(val)
560  class(sll_t_sparse_grid_4d), intent(inout) :: interpolator
561  sll_int32 :: j, l1, l2, l3, level
562  sll_real64 :: val
563  sll_real64, dimension(:), intent(in) :: data, eta
564  sll_real64, dimension(4) :: eta_norm
565  sll_real64, dimension(4) :: phi
566  sll_int32, dimension(4) :: no, l
567  sll_int32, dimension(:, :), allocatable :: ind
568  sll_real64 :: scale
569  sll_int32 :: index
570 
571  sll_allocate(ind(0:interpolator%max_level, 1:4), j)
572 
573  val = 0.0_f64
574  ind(0:1, 1:4) = 0
575 
576  do j = 1, 4
577  eta_norm(j) = (eta(j) - interpolator%eta_min(j))/interpolator%length(j)
578  eta_norm(j) = modulo(eta_norm(j), 1.0_f64)
579 
580  scale = 0.5_f64
581  do level = 2, interpolator%max_level
582  ind(level, j) = ind(level - 1, j)*2
583  if (eta_norm(j) > scale*real(ind(level, j) + 1, f64)) then
584  ind(level, j) = ind(level, j) + 1
585  end if
586  scale = scale*0.5_f64
587  end do
588  end do
589 
590  do level = 0, interpolator%max_level
591  do l1 = 0, min(level, interpolator%levels(1))
592  l(1) = l1
593  no(1) = max(2**(l1 - 1), 1)
594  do l2 = 0, min(level - l1, interpolator%levels(2))
595  l(2) = l2
596  no(2) = max(2**(l2 - 1), 1)
597  do l3 = max(0, level - l1 - l2 - interpolator%levels(4)), min(level - l1 - l2, interpolator%levels(3))
598  no(3) = max(2**(l3 - 1), 1)
599  l(3) = l3
600  l(4) = level - l1 - l2 - l3
601  no(4) = max(2**(l(4) - 1), 1)
602 
603  index = interpolator%index(l1, l2, l3, &
604  l(4)) + ind(l1, 1)*no(2)*no(3)*no(4) &
605  + ind(l2, 2)*no(3)*no(4) + ind(l3, 3)*no(4) + ind(l(4), 4)
606  do j = 1, 4
607  call interpolator%basis_function(real(2**(max(l(j), 1)), f64)*eta_norm(j) &
608  - real(2*ind(l(j), j), f64) - 1.0_f64, phi(j), &
609  interpolator%hierarchy(index)%function_type(j))
610  end do
611  val = val + data(index) &
612  *phi(1)*phi(2)*phi(3)*phi(4)
613  end do
614  end do
615  end do
616  end do
617 
618  end function interpolate_from_hierarchical_surplus
619 
620 !PN DEFINED BUT NOT USED
621 ! function interpolate_from_2D_hierarchical_surplus( interpolator, surplus,eta, dorder,no_in,ind,l ) result(val)
622 ! class(sll_t_sparse_grid_4d), intent(inout) :: interpolator
623 ! sll_real64, dimension(:), intent(in) :: surplus
624 ! sll_int32 :: j,l1,l2,level
625 ! sll_real64 :: val
626 ! sll_real64,dimension(:), intent(in) :: eta
627 ! sll_int32, dimension(:), intent(in) ::dorder
628 ! sll_real64,dimension(4) :: eta_norm
629 ! sll_real64,dimension(4) :: phi
630 ! sll_int32, dimension(:), intent(inout) :: l
631 ! sll_int32, dimension(:), intent(in) :: no_in
632 ! sll_int32, dimension(4) :: no
633 ! sll_int32,dimension(:,:), intent(inout) :: ind
634 ! sll_real64 :: scale
635 ! sll_int32 :: index,maxl2
636 !
637 ! val = 0.0_f64
638 ! ind(1:2,1:4) = 0
639 ! no(dorder(3)) = no_in(dorder(3));
640 ! no(dorder(4)) = no_in(dorder(4));
641 !
642 !! Note this could be organized more efficiently by reusing data
643 ! do j=1,2
644 ! eta_norm(j) = (eta(j)-interpolator%eta_min(dorder(j)))/interpolator%length(dorder(j))
645 ! eta_norm(j) = modulo(eta_norm(j),1.0_f64)
646 ! scale = 0.5_f64
647 ! do level = 0, interpolator%max_level
648 ! ind(level+1,dorder(j)) = ind(level,dorder(j))*2
649 !
650 ! if (eta_norm(j)> scale*(ind(level+1,dorder(j))+1)) then
651 ! ind(level+1,dorder(j)) = ind(level+1,dorder(j))+1
652 ! end if
653 ! scale = scale*0.5_f64
654 ! end do
655 ! end do
656 !
657 !
658 !
659 ! maxl2 = l(dorder(2));
660 ! do l1 = 0, min(interpolator%max_level-l(dorder(3))-l(dorder(4)),interpolator%levels(dorder(1)))
661 ! l(dorder(1)) = l1
662 ! no(dorder(1)) = max(2**(l1-1),1)
663 ! do l2 = 0,min(maxl2,interpolator%max_level-l(dorder(3))-l(dorder(4))-l1)
664 ! l(dorder(2)) = l2
665 ! no(dorder(2)) = max(2**(l2-1),1)
666 ! index = interpolator%index(l(1),l(2),l(3),l(4))+ind(l(1)+1,1)*no(2)*no(3)*no(4)&
667 ! +ind(l(2)+1,2)*no(3)*no(4)+ind(l(3)+1,3)*no(4)+ind(l(4)+1,4)
668 ! do j=1,2
669 ! call interpolator%basis_function(real(2**(max(l(dorder(j)),1)),f64)*eta_norm(j)&
670 ! -real(2*ind(l(dorder(j))+1,dorder(j)),f64)-1.0_f64, phi(j),&
671 ! interpolator%hierarchy(index)%function_type(dorder(j)))
672 ! end do
673 ! val = val + surplus(index)&
674 ! *phi(1)*phi(2)
675 ! !print*, index, val, phi(1)
676 !
677 ! end do
678 ! end do
679 !
680 !
681 ! end function interpolate_from_2D_hierarchical_surplus
682 
683 !!!! End interpolation routines !!!!
684 !------------------------------------------------------------------------------!
685 
686 !------------------------------------------------------------------------------!
687 !!!! Initialization routines !!!!
688 
690  subroutine initialize_sg4d( &
691  interpolator, &
692  levels, &
693  order, &
694  interpolation, &
695  interpolation_type, &
696  eta_min, &
697  eta_max)
698  class(sll_t_sparse_grid_4d), intent(inout) :: interpolator
699  sll_real64, dimension(:), intent(in) :: eta_min
700 
701  sll_real64, dimension(:), intent(in) :: eta_max
702 
703  sll_int32, dimension(:), intent(in) :: levels
704  sll_int32, intent(in) :: order
705  sll_int32, intent(in) :: interpolation
706  sll_int32, intent(in) :: interpolation_type
707  sll_int32 :: i, j, k1, k2, k3, k4, l1, l2, l3, l4, l, counter
708  sll_int32 :: ierr
709  sll_int32, dimension(:), allocatable :: novec, lvec, kvec
710 
711  interpolator%dim = 4;
712  sll_allocate(lvec(interpolator%dim), ierr);
713  sll_allocate(kvec(interpolator%dim), ierr);
714  sll_allocate(novec(interpolator%dim), ierr);
715  interpolator%max_level = levels(1);
716  do l = 2, interpolator%dim
717  interpolator%max_level = max(levels(l), interpolator%max_level);
718  end do
719 ! if (interpolator%modified == 1) then
720 ! interpolator%max_level = interpolator%max_level + 2;
721 ! end if
722 
723  interpolator%size_basis = 0;
724  do l = 0, interpolator%max_level
725  do l1 = 0, min(l, levels(1))
726  do l2 = 0, min(l - l1, levels(2))
727  do l3 = max(0, l - l1 - l2 - levels(4)), min(l - l1 - l2, levels(3))
728  l4 = l - l1 - l2 - l3
729  interpolator%size_basis = &
730  interpolator%size_basis + &
731  max(2**(l1 - 1), 1)*max(2**(l2 - 1), 1)*max(2**(l3 - 1), 1)*max(2**(l4 - 1), 1)
732  end do
733  end do
734  end do
735  end do
736 
737  call interpolator%initialize_sg(levels, order, interpolation, &
738  interpolation_type, eta_min, eta_max);
739  sll_allocate(interpolator%index(0:interpolator%levels(1),0:interpolator%levels(2),0:interpolator%levels(3),0:interpolator%levels(4)),ierr)
740 
741  ! Set the hierarchy of the grid
742  counter = 1
743  do l = 0, interpolator%max_level
744  interpolator%level_mapping(l) = counter;
745  do l1 = 0, min(l, interpolator%levels(1))
746  novec(1) = max(2**(l1 - 1), 1)
747  lvec(1) = l1;
748  do l2 = 0, min(l - l1, interpolator%levels(2))
749  novec(2) = max(2**(l2 - 1), 1)
750  lvec(2) = l2;
751  do l3 = max(0, l - l1 - l2 - interpolator%levels(4)), min(l - l1 - l2, interpolator%levels(3))
752  novec(3) = max(2**(l3 - 1), 1)
753  lvec(3) = l3
754  lvec(4) = l - l1 - l2 - l3
755  l4 = lvec(4)
756  novec(4) = max(2**(l4 - 1), 1)
757  interpolator%index(l1, l2, l3, l4) = counter
758  do k1 = 0, novec(1) - 1
759  kvec(1) = k1
760  do k2 = 0, novec(2) - 1
761  kvec(2) = k2
762  do k3 = 0, novec(3) - 1
763  kvec(3) = k3
764  do k4 = 0, novec(4) - 1
765  kvec(4) = k4;
766  do j = 1, interpolator%dim
767  call set_hierarchy_info(interpolator, counter, j, &
768  lvec, kvec, novec);
769  end do
770  counter = counter + 1;
771  end do
772  end do
773  end do
774  end do
775  end do
776  end do
777  end do
778  end do
779  interpolator%level_mapping(interpolator%max_level + 1) = counter;
780  ! Now rescale all the coordinates to the actual mesh size
781  do i = 1, interpolator%size_basis
782  do j = 1, interpolator%dim
783  interpolator%hierarchy(i)%coordinate(j) = &
784  interpolator%hierarchy(i)%coordinate(j)* &
785  interpolator%length(j) + &
786  interpolator%eta_min(j)
787  end do
788  end do
789 
790  end subroutine initialize_sg4d
791 
793  subroutine set_hierarchy_info(interpolator, counter, cdim, lvecin, kvecin, novecin)
794  class(sll_t_sparse_grid_4d), intent(inout) :: interpolator
795  sll_int32 :: ld
796  sll_int32 :: kd
797  sll_int32, intent(in) :: cdim
798  sll_int32, intent(in) :: counter
799  sll_int32, dimension(:), intent(in) :: lvecin
800  sll_int32, dimension(:), intent(in) :: kvecin
801  sll_int32, dimension(:), intent(in) :: novecin
802 
803  sll_int32, dimension(4) :: lvec, kvec, novec
804  sll_int32 :: jj, stride
805 
806  do jj = 1, interpolator%dim
807  lvec(jj) = lvecin(jj);
808  kvec(jj) = kvecin(jj);
809  novec(jj) = novecin(jj);
810  end do
811  ld = lvec(cdim);
812  kd = kvec(cdim);
813  interpolator%hierarchy(counter)%level(cdim) = ld;
814  interpolator%hierarchy(counter)%index_on_level(cdim) = kd;
815  stride = cdim*2 - 1
816  if (ld == 0) then
817  interpolator%hierarchy(counter)%coordinate(cdim) = 0.0_f64
818  interpolator%hierarchy(counter)%parent(stride) = &
819  counter
820  interpolator%hierarchy(counter)%parent(stride + 1) = &
821  counter
822  interpolator%hierarchy(counter)%function_type(cdim) = 0
823  else
824  interpolator%hierarchy(counter)%coordinate(cdim) = &
825  1.0_f64/(2.0_f64**ld) + real(kd, f64)*1.0_f64/(2.0_f64**(ld - 1))
826 
827  lvec(cdim) = lvec(cdim) - 1;
828  novec(cdim) = max(novec(cdim)/2, 1);
829  ! This one is actually only a neighboring point not the direct parent.
830  kvec(cdim) = modulo((kd + 1)/2, max(2**(ld - 2), 1));
831  interpolator%hierarchy(counter)%parent( &
832  modulo(kd, 2) + stride) = &
833  interpolator%hierarchy( &
834  interpolator%index(lvec(1), lvec(2), lvec(3), lvec(4)) + &
835  kvec(1)*novec(2)*novec(3)*novec(4) + &
836  kvec(2)*novec(3)*novec(4) + &
837  kvec(3)*novec(4) + kvec(4))%parent(stride)
838  ! This is the actual parent.
839  kvec(cdim) = kd/2;
840  interpolator%hierarchy(counter)%parent( &
841  modulo(kd + 1, 2) + stride) = &
842  interpolator%index(lvec(1), lvec(2), lvec(3), lvec(4)) + &
843  kvec(1)*novec(2)*novec(3)*novec(4) + kvec(2)*novec(3)*novec(4) + &
844  kvec(3)*novec(4) + kvec(4)
845  ! Now tell my parent that I am his child.
846  interpolator%hierarchy( &
847  interpolator%hierarchy(counter)%parent( &
848  modulo(kd + 1, 2) + stride))%children(modulo(kd, 2) + stride) = counter
849  if (ld == 1) then
850  interpolator%hierarchy( &
851  interpolator%hierarchy(counter)%parent( &
852  modulo(kd + 1, 2) + stride))%children(modulo(kd + 1, 2) + stride) = counter
853  end if
854  if (interpolator%order == 1) then
855  interpolator%hierarchy(counter)% &
856  function_type(cdim) = -1
857  elseif (ld == 1 .OR. interpolator%order == 2) then
858  interpolator%hierarchy(counter)% &
859  function_type(cdim) = 1 + modulo(kd, 2)
860  elseif (ld == 2 .OR. interpolator%order == 3) then !(order==3) then!
861  interpolator%hierarchy(counter)% &
862  function_type(cdim) = 3 + modulo(kd, 4)
863  else
864  interpolator%hierarchy(counter)% &
865  function_type(cdim) = 7 + modulo(kd, 8)
866  end if
867  end if
868 
869  end subroutine set_hierarchy_info
870 
871 !!!! End initialization routines !!!!
872 !------------------------------------------------------------------------------!
873 
874 !------------------------------------------------------------------------------!
875 !PN DEFINED BUT NOT USED
876 !!> Functions to evaluate fg on sg and sg on fg
877 !subroutine fg_to_sg(interpolator,fg_values,sg_values)
878 !sll_real64, dimension(:,:,:,:), intent(in) :: fg_values
879 !sll_real64, dimension(:), intent(out) :: sg_values
880 !class(sll_t_sparse_grid_4d), intent(in) :: interpolator
881 !sll_int32 :: j
882 !sll_int32, dimension(4) :: fg_ind
883 !
884 !do j=1,interpolator%size_basis
885 ! fg_ind = fg_index(interpolator,j);
886 ! sg_values(j) = fg_values(fg_ind(1),fg_ind(2),fg_ind(3), fg_ind(4));
887 !end do
888 !
889 !end subroutine fg_to_sg
890 
892  function fg_index(interpolator, sg_index)
893  sll_int32, intent(in) :: sg_index
894  sll_int32, dimension(4) :: fg_index
895  class(sll_t_sparse_grid_4d), intent(in) :: interpolator
896  sll_int32 :: j
897 
898  do j = 1, interpolator%dim
899  fg_index(j) = 2**(interpolator%levels(j) - &
900  interpolator%hierarchy(sg_index)%level(j))* &
901  (1 + 2*interpolator%hierarchy(sg_index)%index_on_level(j)) + 1;
902  end do
903 
904  end function fg_index
905 
906 ! End functions fg_to_sg and sg_to_fg
907 !------------------------------------------------------------------------------!
908 
909 end module sll_m_sparse_grid_4d
sll_m_constants
Fortran module where set some physical and mathematical constants.
Definition: sll_m_constants.F90:23
sll_m_constants::sll_p_pi
real(kind=f64), parameter, public sll_p_pi
Definition: sll_m_constants.F90:51
sll_m_sparse_grid_4d
Implementation of a 4D sparse grid with interpolation routines.
Definition: sll_m_sparse_grid_4d.F90:6
sll_m_sparse_grid_4d::interpolate_disp_nconst_in_1d
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....
Definition: sll_m_sparse_grid_4d.F90:131
sll_m_sparse_grid_4d::initialize_sg4d
subroutine initialize_sg4d(interpolator, levels, order, interpolation, interpolation_type, eta_min, eta_max)
Initialization function. Set up the hierarchy of the sparse grid.
Definition: sll_m_sparse_grid_4d.F90:698
sll_m_sparse_grid_4d::interpolate_const_disp
subroutine interpolate_const_disp(interpolator, dorder, displacement, data_in, data_out, hiera)
Interpolation function for interpolation at (constantly) displaced grid points; displacement only in ...
Definition: sll_m_sparse_grid_4d.F90:61
sll_m_sparse_grid_4d::fg_index
integer(kind=i32) function, dimension(4) 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...
Definition: sll_m_sparse_grid_4d.F90:893
sll_m_sparse_grid_4d::interpolate_from_interpolant_value
real(kind=f64) function interpolate_from_interpolant_value(interpolator, data, eta)
Compute the value of the sparse grid interpolant at position eta.
Definition: sll_m_sparse_grid_4d.F90:46
sll_m_sparse_grid_4d::interpolate4d_disp_linnconst_in_1d
subroutine interpolate4d_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))
Definition: sll_m_sparse_grid_4d.F90:232
sll_m_sparse_grid_4d::set_hierarchy_info
subroutine set_hierarchy_info(interpolator, counter, cdim, lvecin, kvecin, novecin)
For a given sparse grid point fill the hierarchy information (4D specific)
Definition: sll_m_sparse_grid_4d.F90:794
sll_m_sparse_grid_4d::interpolate_disp_nconst_in_2d
subroutine interpolate_disp_nconst_in_2d(interpolator, displacement, dorder, data_in, data_out)
As previous function but with displacement displacement*coordinate(dorder(2))
Definition: sll_m_sparse_grid_4d.F90:332
sll_m_sparse_grid_4d::interpolate_from_hierarchical_surplus
real(kind=f64) function interpolate_from_hierarchical_surplus(interpolator, data, eta)
Definition: sll_m_sparse_grid_4d.F90:560
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_4d::sll_t_sparse_grid_4d
Sparse grid object for 4d with interpolation routines. Note in 4d we have only an implementation of a...
Definition: sll_m_sparse_grid_4d.F90:26
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
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:39 for Semi-Lagrangian Library by   doxygen 1.9.1