Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_splines_pp.F90
Go to the documentation of this file.
1 
7 
8 module sll_m_splines_pp
9 
10  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
11 #include "sll_assert.h"
12 #include "sll_errors.h"
13 #include "sll_memory.h"
14 #include "sll_working_precision.h"
15 
16  implicit none
17 
18  public :: &
19  sll_t_spline_pp_1d, &
20  sll_t_spline_pp_2d, &
21  sll_t_spline_pp_3d, &
22  sll_s_spline_pp_b_to_pp_1d, &
23  sll_s_spline_pp_b_to_pp_2d, &
24  sll_s_spline_pp_b_to_pp_2d_periodic, &
25  sll_s_spline_pp_b_to_pp_3d, &
26  sll_s_spline_pp_b_to_pp_3d_clamped, &
27  sll_s_spline_pp_b_to_pp_3d_clamped_2full, &
28  sll_s_spline_pp_horner_m_1d, &
29  sll_s_spline_pp_horner_m_2d, &
30  sll_s_spline_pp_horner_m_3d, &
31  sll_s_spline_pp_horner_primitive_1d, &
32  sll_f_spline_pp_horner_1d, &
33  sll_f_spline_pp_horner_2d, &
34  sll_f_spline_pp_horner_3d, &
35  sll_f_spline_pp_horner_3d_d1, &
36  sll_f_spline_pp_horner_3d_d2, &
37  sll_f_spline_pp_horner_3d_d3, &
38  sll_f_spline_pp_horner_derivative_3d, &
39  sll_s_spline_pp_init_1d, &
40  sll_s_spline_pp_init_2d, &
41  sll_s_spline_pp_init_3d, &
42  sll_s_spline_pp_free_1d, &
43  sll_s_spline_pp_free_2d, &
44  sll_s_spline_pp_free_3d, &
45  sll_s_spline_pp_pp_to_b_1d, &
46  sll_s_spline_evaluate_basis_b_form_1d_clamped, &
47  sll_s_spline_evaluate_basis_b_form_1d_periodic, &
48  sll_s_spline_evaluate_basis_b_form_3d_clamped, &
49  sll_s_spline_evaluate_basis_b_form_3d_periodic, &
50  sll_p_boundary_periodic, &
51  sll_p_boundary_clamped, &
52  sll_p_boundary_clamped_square, &
53  sll_p_boundary_clamped_cubic, &
54  sll_p_boundary_clamped_clampeddiri, &
55  sll_p_boundary_clampeddiri, &
56  sll_p_boundary_clampeddiri_clamped
57 
58  private
59  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
61  sll_int32, parameter :: sll_p_boundary_periodic = 0
62  sll_int32, parameter :: sll_p_boundary_clamped = 1
63  sll_int32, parameter :: sll_p_boundary_clamped_clampeddiri = 2
64  sll_int32, parameter :: sll_p_boundary_clampeddiri = 3
65  sll_int32, parameter :: sll_p_boundary_clampeddiri_clamped = 4
66  sll_int32, parameter :: sll_p_boundary_clamped_square = 5
67  sll_int32, parameter :: sll_p_boundary_clamped_cubic = 6
68 
69  sll_real64, parameter :: inv_2 = 1._f64/2._f64
70  sll_real64, parameter :: inv_3 = 1._f64/3._f64
71  sll_real64, parameter :: inv_4 = 1._f64/4._f64
72  sll_real64, parameter :: inv_6 = 1._f64/6._f64
73  sll_real64, parameter :: inv_7 = 1._f64/7._f64
74  sll_real64, parameter :: inv_8 = 1._f64/8._f64
75  sll_real64, parameter :: inv_12 = 1._f64/12._f64
76  sll_real64, parameter :: inv_16 = 1._f64/16._f64
77  sll_real64, parameter :: inv_18 = 1._f64/18._f64
78  sll_real64, parameter :: inv_20 = 1._f64/20._f64
79  sll_real64, parameter :: inv_24 = 1._f64/24._f64
80  sll_real64, parameter :: inv_30 = 1._f64/30._f64
81  sll_real64, parameter :: inv_36 = 1._f64/36._f64
82  sll_real64, parameter :: inv_48 = 1._f64/48._f64
83  sll_real64, parameter :: inv_60 = 1._f64/60._f64
84  sll_real64, parameter :: inv_72 = 1._f64/72._f64
85  sll_real64, parameter :: inv_80 = 1._f64/80._f64
86  sll_real64, parameter :: inv_120 = 1._f64/120._f64
87  sll_real64, parameter :: inv_240 = 1._f64/240._f64
88  sll_real64, parameter :: inv_144 = 1._f64/144._f64
89  ! sll_real64, parameter :: inv_600 = 1._f64/600._f64
90  sll_real64, parameter :: inv_720 = 1._f64/720._f64
91  sll_real64, parameter :: inv_5040 = 1._f64/5040._f64
93  type :: sll_t_spline_pp_1d
94 
95  sll_int32 :: degree
96  sll_real64, allocatable :: poly_coeffs(:, :)
97  sll_real64, allocatable :: poly_coeffs_fp(:, :)
98  sll_real64, allocatable :: poly_coeffs_fpa(:, :)
99  sll_real64, allocatable :: poly_coeffs_fd(:, :)
100  sll_int32 :: n_cells
101  sll_int32 :: n_coeffs
102  sll_real64, allocatable :: scratch_b(:)
103  sll_real64, allocatable :: scratch_p(:)
104 
105  sll_real64, allocatable :: poly_coeffs_boundary_left(:, :, :)
106  sll_real64, allocatable :: poly_coeffs_boundary_right(:, :, :)
107  sll_real64, allocatable :: poly_coeffs_fd_boundary_left(:, :, :)
108  sll_real64, allocatable :: poly_coeffs_fd_boundary_right(:, :, :)
109  sll_int32 :: boundary_conditions = sll_p_boundary_periodic
110 
111  end type sll_t_spline_pp_1d
113  type :: sll_t_spline_pp_2d
114  type(sll_t_spline_pp_1d) spline1
115  type(sll_t_spline_pp_1d) spline2
116 
117  end type sll_t_spline_pp_2d
118 
120  type :: sll_t_spline_pp_3d
121  type(sll_t_spline_pp_1d) spline1
122  type(sll_t_spline_pp_1d) spline2
123  type(sll_t_spline_pp_1d) spline3
124 
125  sll_real64, allocatable :: scratch_b(:)
126  sll_real64, allocatable :: scratch_p(:)
127  sll_real64, allocatable :: scratch_pp(:)
128  sll_real64, allocatable :: scratch_coeffs(:, :)
129 
130  end type sll_t_spline_pp_3d
131 
132 contains
133 
135  subroutine sll_s_spline_pp_pp_to_b_1d(self, n_cells, pp_coeffs, b_coeffs)
136  type(sll_t_spline_pp_1d), intent(in):: self
137  sll_int32, intent(in) :: n_cells
138  sll_real64, intent(out) :: b_coeffs(0:n_cells - 1)
139  sll_real64, intent(in):: pp_coeffs(self%degree + 1, 0:n_cells - 1)
140 
141  sll_int32 :: i, j, imin
142 
143  b_coeffs = 0.0_f64
144  do i = 0, n_cells - 1
145  imin = i - self%degree
146  do j = 0, self%degree
147  b_coeffs(modulo(imin + j, n_cells)) = b_coeffs(modulo(imin + j, n_cells)) + &
148  sum(pp_coeffs(:, i)*self%poly_coeffs(:, j + 1))
149  end do
150  end do
151 
152  end subroutine sll_s_spline_pp_pp_to_b_1d
153 
155  subroutine sll_s_spline_pp_b_to_pp_1d_clamped(self, n_cells, b_coeffs, pp_coeffs)
156  type(sll_t_spline_pp_1d), intent(in):: self
157  sll_int32, intent(in) :: n_cells
158  sll_real64, intent(in) :: b_coeffs(n_cells + self%degree)
159  sll_real64, intent(out):: pp_coeffs(self%degree + 1, n_cells)
160  sll_int32 :: i, upper
161 
162  do i = 1, self%degree - 1
163  call sll_s_spline_pp_b_to_pp_1d_cella(self%degree, self%poly_coeffs_boundary_left(:, :, i), &
164  b_coeffs(i:i + self%degree), pp_coeffs(:, i))
165  end do
166 
167  upper = n_cells - self%degree + 1
168  do i = self%degree, upper
169  call sll_s_spline_pp_b_to_pp_1d_cella(self%degree, self%poly_coeffs, &
170  b_coeffs(i:i + self%degree), pp_coeffs(:, i))
171  end do
172 
173  do i = upper + 1, n_cells
174  call sll_s_spline_pp_b_to_pp_1d_cella(self%degree, self%poly_coeffs_boundary_right(:, :, i - upper), &
175  b_coeffs(i:i + self%degree), pp_coeffs(:, i))
176  end do
177 
178  end subroutine sll_s_spline_pp_b_to_pp_1d_clamped
179 
181  subroutine sll_s_spline_pp_b_to_pp_1d_clamped_clampeddiri(self, n_cells, b_coeffs, pp_coeffs)
182  type(sll_t_spline_pp_1d), intent(in):: self
183  sll_int32, intent(in) :: n_cells
184  sll_real64, intent(in) :: b_coeffs(n_cells + self%degree - 1)
185  sll_real64, intent(out):: pp_coeffs(self%degree + 1, n_cells)
186  sll_int32 :: i, upper
187 
188  do i = 1, self%degree - 1
189  call sll_s_spline_pp_b_to_pp_1d_cella(self%degree, self%poly_coeffs_boundary_left(:, :, i), &
190  b_coeffs(i:i + self%degree), pp_coeffs(:, i))
191  end do
192 
193  upper = n_cells - self%degree + 1
194  do i = self%degree, upper
195  call sll_s_spline_pp_b_to_pp_1d_cella(self%degree, self%poly_coeffs, &
196  b_coeffs(i:i + self%degree), pp_coeffs(:, i))
197  end do
198 
199  do i = upper + 1, n_cells - 1
200  call sll_s_spline_pp_b_to_pp_1d_cella(self%degree, self%poly_coeffs_boundary_right(:, :, i - upper), &
201  b_coeffs(i:i + self%degree), pp_coeffs(:, i))
202  end do
203  call sll_s_spline_pp_b_to_pp_1d_cella(self%degree, self%poly_coeffs_boundary_right(:, :, i - upper), &
204  [b_coeffs(i:i + self%degree - 1), 0.0_f64], pp_coeffs(:, i))
205 
206  end subroutine sll_s_spline_pp_b_to_pp_1d_clamped_clampeddiri
207 
209  subroutine sll_s_spline_pp_b_to_pp_1d(self, n_cells, b_coeffs, pp_coeffs)
210  type(sll_t_spline_pp_1d), intent(in):: self
211  sll_int32, intent(in) :: n_cells
212  sll_real64, intent(in) :: b_coeffs(n_cells)
213  sll_real64, intent(out):: pp_coeffs(self%degree + 1, n_cells)
214  sll_int32 :: i
215 
216  do i = 1, self%degree
217  call sll_s_spline_pp_b_to_pp_1d_cell(self, (/b_coeffs(n_cells - self%degree + i:n_cells), b_coeffs(1:i)/), pp_coeffs(:, i))
218  end do
219 
220  do i = self%degree + 1, n_cells
221  call sll_s_spline_pp_b_to_pp_1d_cell(self, b_coeffs(i - self%degree:i), pp_coeffs(:, i))
222  end do
223 
224  end subroutine sll_s_spline_pp_b_to_pp_1d
225 
227  subroutine sll_s_spline_pp_b_to_pp_1d_cell(self, b_coeffs, pp_coeffs)
228  type(sll_t_spline_pp_1d), intent(in):: self
229  sll_real64, intent(in) :: b_coeffs(self%degree + 1)
230  sll_real64, intent(out) :: pp_coeffs(self%degree + 1)
231  sll_int32 :: i
232  sll_int32 :: j
233  sll_int32 :: degp1
234  degp1 = self%degree + 1
235  pp_coeffs = 0.0_f64
236  do i = 1, degp1
237  do j = 1, degp1
238  pp_coeffs(j) = pp_coeffs(j) + b_coeffs(i)*self%poly_coeffs(j, i)
239  end do
240  end do
241 
242  end subroutine sll_s_spline_pp_b_to_pp_1d_cell
243 
245  subroutine sll_s_spline_pp_b_to_pp_1d_cella(degree, pp_b, b_coeffs, pp_coeffs)
246  sll_int32, intent(in) :: degree
247  sll_real64, intent(in) :: pp_b(degree + 1, degree + 1)
248  sll_real64, intent(in) :: b_coeffs(degree + 1)
249  sll_real64, intent(out) :: pp_coeffs(degree + 1)
250  sll_int32 :: i
251  sll_int32 :: j
252  sll_int32 :: degp1
253  degp1 = degree + 1
254  pp_coeffs = 0.0_f64
255  do i = 1, degp1
256  do j = 1, degp1
257  pp_coeffs(j) = pp_coeffs(j) + b_coeffs(i)*pp_b(j, i)
258  end do
259  end do
260 
261  end subroutine sll_s_spline_pp_b_to_pp_1d_cella
262 
265  subroutine sll_s_spline_pp_b_to_pp_2d_periodic(self, n_cells, b_coeffs, pp_coeffs)
266  type(sll_t_spline_pp_2d), intent(inout):: self
267  sll_int32, intent(in) :: n_cells(2)
268  sll_real64, intent(in) :: b_coeffs(n_cells(1)*n_cells(2))
269  sll_real64, intent(out):: pp_coeffs((self%spline1%degree + 1)*(self%spline2%degree + 1), n_cells(1)*n_cells(2))
270  sll_int32 :: i, j
271  sll_int32 :: degree1, degree2
272  degree1 = self%spline1%degree
273  degree2 = self%spline2%degree
274 
275  do j = 1, n_cells(2)
276  do i = 1, n_cells(1)
277  call sll_s_spline_pp_b_to_pp_2d_cell(self%spline1, self%spline2, n_cells, b_coeffs, pp_coeffs, i, j)
278  end do
279  end do
280 
281  end subroutine sll_s_spline_pp_b_to_pp_2d_periodic
282 
284  subroutine sll_s_spline_pp_b_to_pp_2d(self, n_cells, b_coeffs, pp_coeffs)
285  type(sll_t_spline_pp_2d), intent(inout):: self
286  sll_int32, intent(in) :: n_cells(2)
287  sll_real64, intent(in) :: b_coeffs(:)!(n_cells(1)*n_cells(2)) !< coefficients of spline in B-form
288  sll_real64, intent(out):: pp_coeffs((self%spline1%degree + 1)*(self%spline2%degree + 1), n_cells(1)*n_cells(2))
289  sll_int32 :: i, j, l1, l2
290  sll_int32 :: degree1, degree2
291  sll_real64 :: pp_coeffs_local(self%spline1%degree + 1, self%spline2%degree + 1)
292  sll_int32 :: offset1, offset2, cell, n_coeffs(2), index, upper(2)
293 
294  degree1 = self%spline1%degree
295  degree2 = self%spline2%degree
296  n_coeffs(1) = self%spline1%n_coeffs
297  n_coeffs(2) = self%spline2%n_coeffs
298 
299  select case (self%spline1%boundary_conditions)
300  case (sll_p_boundary_periodic)
301  offset1 = -degree1
302  case (sll_p_boundary_clampeddiri:sll_p_boundary_clamped_cubic)
303  offset1 = -1
304  case default
305  offset1 = 0
306  end select
307 
308  upper(1) = n_cells(1) - degree1 + 1
309  upper(2) = n_cells(2) - degree2 + 1
310  do j = 1, n_cells(2)
311  do i = 1, n_cells(1)
312  cell = (j - 1)*n_cells(1) + i
313  do l2 = 0, degree2
314  select case (self%spline2%boundary_conditions)
315  case (sll_p_boundary_periodic)
316  offset2 = modulo(j - degree2 + l2 - 1, n_coeffs(2))*n_coeffs(1) ! periodic
317  case (sll_p_boundary_clampeddiri)
318  offset2 = (j + l2 - 2)*n_coeffs(1)
319  if (j == 1 .and. l2 == 0) then
320  pp_coeffs_local(:, 1) = 0.0_f64
321  !print*, 'a'
322  !continue
323  !goto 100
324  cycle
325  elseif (j == n_cells(2) .and. l2 == degree2) then
326  pp_coeffs_local(:, degree2 + 1) = 0.0_f64
327  !continue
328  !goto 100
329  cycle
330  end if
331  case (sll_p_boundary_clamped_clampeddiri)
332  offset2 = (j + l2 - 1)*n_coeffs(1)
333  if (j == n_cells(2) .and. l2 == degree2) then
334  pp_coeffs_local(:, degree2 + 1) = 0.0_f64
335  !continue
336  !goto 100
337  cycle
338  end if
339  case (sll_p_boundary_clampeddiri_clamped)
340  offset2 = (j + l2 - 2)*n_coeffs(1)
341  if (j == 1 .and. l2 == 0) then
342  pp_coeffs_local(:, 1) = 0.0_f64
343  !continue
344  !goto 100
345  cycle
346  end if
347  !case ( sll_p_boundary_clampeddiri:sll_p_boundary_clampeddiri_clamped)
348  ! offset2 = (j+l2-2)*n_coeffs(1)
349  case default
350  offset2 = (j + l2 - 1)*n_coeffs(1)
351  end select
352 
353  do l1 = 0, degree1
354  select case (self%spline1%boundary_conditions)
355  case (sll_p_boundary_periodic)
356  index = modulo(i + offset1 + l1 - 1, n_coeffs(1)) + 1
357  !print*, i, j, l1, l2, offset2, index
358  self%spline1%scratch_b(l1 + 1) = b_coeffs(offset2 + index)
359  case (sll_p_boundary_clamped)
360  index = i + offset1 + l1
361  self%spline1%scratch_b(l1 + 1) = b_coeffs(offset2 + index)
362  ! Set to zero for Dirichlet if we are at the boundary
363  case (sll_p_boundary_clamped_square)
364  index = i + offset1 + l1
365  self%spline1%scratch_b(l1 + 1) = b_coeffs(offset2 + index)
366  ! Set to zero for Dirichlet if we are at the boundary
367  case (sll_p_boundary_clamped_cubic)
368  index = i + offset1 + l1
369  self%spline1%scratch_b(l1 + 1) = b_coeffs(offset2 + index)
370  ! Set to zero for Dirichlet if we are at the boundary
371  case (sll_p_boundary_clamped_clampeddiri)
372  if (i == n_cells(1) .and. l1 == degree1) then
373  self%spline1%scratch_b(l1 + 1) = 0.0_f64
374  else
375  index = i + offset1 + l1
376  self%spline1%scratch_b(l1 + 1) = b_coeffs(offset2 + index)
377  end if
378  case (sll_p_boundary_clampeddiri_clamped)
379  if (i == 1 .and. l1 == 0) then
380  self%spline1%scratch_b(l1 + 1) = 0.0_f64
381  else
382  index = i + offset1 + l1
383  self%spline1%scratch_b(l1 + 1) = b_coeffs(offset2 + index)
384  end if
385  case (sll_p_boundary_clampeddiri)
386  if (i == 1 .and. l1 == 0) then
387  self%spline1%scratch_b(l1 + 1) = 0.0_f64
388  elseif (i == n_cells(1) .and. l1 == degree1) then
389  self%spline1%scratch_b(l1 + 1) = 0.0_f64
390  else
391  index = i + offset1 + l1
392  self%spline1%scratch_b(l1 + 1) = b_coeffs(offset2 + index)
393  end if
394  end select
395  end do
396  ! For clamped splines, we need to use the boundary coefficients
397  select case (self%spline1%boundary_conditions)
398  case (sll_p_boundary_periodic)
399  call sll_s_spline_pp_b_to_pp_1d_cella(degree1, self%spline1%poly_coeffs, &
400  self%spline1%scratch_b, pp_coeffs_local(:, l2 + 1))
401  case default
402  if (i > degree1 - 1) then
403  if (i < upper(1) + 1) then
404  call sll_s_spline_pp_b_to_pp_1d_cella(degree1, self%spline1%poly_coeffs, &
405  self%spline1%scratch_b, pp_coeffs_local(:, l2 + 1))
406  else
407  call sll_s_spline_pp_b_to_pp_1d_cella(degree1, &
408  self%spline1%poly_coeffs_boundary_right(:, :, i - upper(1)), &
409  self%spline1%scratch_b, pp_coeffs_local(:, l2 + 1))
410  end if
411  else
412 
413  call sll_s_spline_pp_b_to_pp_1d_cella(degree1, &
414  self%spline1%poly_coeffs_boundary_left(:, :, i), &
415  self%spline1%scratch_b, pp_coeffs_local(:, l2 + 1))
416  end if
417 
418  end select
419  !call sll_s_spline_pp_b_to_pp_1d_cell(self%spline1, &
420  ! self%spline1%scratch_b,pp_coeffs_local(:,l2+1))
421  !100 continue
422  end do
423  do l1 = 0, degree1
424  self%spline2%scratch_b = pp_coeffs_local(l1 + 1, :)
425  select case (self%spline2%boundary_conditions)
426  case (sll_p_boundary_periodic)
427  call sll_s_spline_pp_b_to_pp_1d_cella(degree2, self%spline2%poly_coeffs, &
428  self%spline2%scratch_b, self%spline2%scratch_p)
429  case default
430  if (j > degree2 - 1) then
431  if (j < upper(2) + 1) then
432  call sll_s_spline_pp_b_to_pp_1d_cella(degree2, &
433  self%spline2%poly_coeffs, &
434  self%spline2%scratch_b, self%spline2%scratch_p)
435  else
436  call sll_s_spline_pp_b_to_pp_1d_cella(degree2, &
437  self%spline2%poly_coeffs_boundary_right(:, :, j - upper(2)), &
438  self%spline2%scratch_b, self%spline2%scratch_p)
439  end if
440  else
441 
442  call sll_s_spline_pp_b_to_pp_1d_cella(degree2, &
443  self%spline2%poly_coeffs_boundary_left(:, :, j), &
444  self%spline2%scratch_b, self%spline2%scratch_p)
445  end if
446 
447  end select
448 
449  !call sll_s_spline_pp_b_to_pp_1d_cell(self%spline2, &
450  ! self%spline2%scratch_b,self%spline2%scratch_p)
451  do l2 = 0, degree2
452  pp_coeffs(l2*(degree1 + 1) + l1 + 1, cell) = self%spline2%scratch_p(l2 + 1)
453  end do
454  end do
455  end do
456  end do
457 
458  end subroutine sll_s_spline_pp_b_to_pp_2d
459 
460 !!$ !> Convert 2d spline in B form in a cell to spline in pp form with periodic boundary conditions
461 !!$ subroutine sll_s_spline_pp_b_to_pp_2d_cell_clamped(spline1,spline2,n_cells, b_coeffs, pp_coeffs,i,j)
462 !!$ type( sll_t_spline_pp_1d), intent(inout):: spline1 !< arbitrary degree 1d spline
463 !!$ type( sll_t_spline_pp_1d), intent(inout):: spline2 !< arbitrary degree 1d spline
464 !!$ sll_int32, intent(in) :: n_cells(2) !< number of gridcells
465 !!$ sll_real64,intent(in) :: b_coeffs(n_cells(1)*n_cells(2)) !< coefficients of spline in B-form
466 !!$ sll_real64,intent(inout) :: pp_coeffs((spline1%degree+1)*(spline2%degree+1),n_cells(1)*n_cells(2)) !< coefficients of spline in pp-form
467 !!$
468 !!$ !> convert b-coefficients to pp-coefficients in first dimension
469 !!$ do j=1, n_cells(2)
470 !!$ do i=1, n_cells(1)
471 !!$ do l=0, degree(2)
472 !!$ spline1%scratch_b=b_coeffs(i-degree1+(j-degp2+l)*n_cells(1):i+(j-degp2+l)*n_cells(1))
473 !!$
474 !!$ end subroutine sll_s_spline_pp_b_to_pp_2d_cell_clamped
475 
477  subroutine sll_s_spline_pp_b_to_pp_2d_cell(spline1, spline2, n_cells, b_coeffs, pp_coeffs, i, j)
478  type(sll_t_spline_pp_1d), intent(inout):: spline1
479  type(sll_t_spline_pp_1d), intent(inout):: spline2
480  sll_int32, intent(in) :: n_cells(2)
481  sll_real64, intent(in) :: b_coeffs(n_cells(1)*n_cells(2))
482  sll_real64, intent(inout) :: pp_coeffs((spline1%degree + 1)*(spline2%degree + 1), n_cells(1)*n_cells(2))
483  sll_int32, intent(in) :: i, j
484  sll_int32 :: k, l
485  sll_int32 :: degree1, degree2, degp1, degp2
486  degree1 = spline1%degree
487  degree2 = spline2%degree
488  degp1 = degree1 + 1
489  degp2 = degree2 + 1
491  if (i > degree1) then
492  if (j > degree2) then
493  do l = 0, degree2
494  spline1%scratch_b = b_coeffs(i - degree1 + (j - degp2 + l)*n_cells(1):i + (j - degp2 + l)*n_cells(1))
495  call sll_s_spline_pp_b_to_pp_1d_cell(spline1, spline1%scratch_b, pp_coeffs(1 + l*degp1:degp1*(l + 1), i + n_cells(1)*(j - 1)))
496  end do
497  else
499  do l = 0, degree2
500  spline1%scratch_b = b_coeffs(i - degree1 + modulo(j - degp2 + l, n_cells(2))*n_cells(1):i + modulo(j - degp2 + l, n_cells(2))*n_cells(1))
501  call sll_s_spline_pp_b_to_pp_1d_cell(spline1, spline1%scratch_b, pp_coeffs(1 + l*degp1:degp1*(l + 1), i + n_cells(1)*(j - 1)))
502  end do
503  end if
504  else
506  do l = 0, degree2
507  do k = 0, degree1
508  spline1%scratch_b(k + 1) = b_coeffs(modulo(i - degp1 + k, n_cells(1)) + 1 + modulo(j - degp2 + l, n_cells(2))*n_cells(1))
509  end do
510 
511  call sll_s_spline_pp_b_to_pp_1d_cell(spline1, spline1%scratch_b, pp_coeffs(1 + l*degp1:degp1*(l + 1), i + n_cells(1)*(j - 1)))
512  end do
513  end if
514 
516  do k = 1, degp1
517  do l = 1, degp2
518  spline2%scratch_b(l) = pp_coeffs(k + degp1*(l - 1), i + n_cells(1)*(j - 1))
519  end do
520  call sll_s_spline_pp_b_to_pp_1d_cell(spline2, spline2%scratch_b, spline2%scratch_p)
521  do l = 1, degp2
522  pp_coeffs(k + degp1*(l - 1), i + n_cells(1)*(j - 1)) = spline2%scratch_p(l)
523  end do
524  end do
525  end subroutine sll_s_spline_pp_b_to_pp_2d_cell
526 
528  subroutine sll_s_spline_pp_b_to_pp_3d(self, n_cells, b_coeffs, pp_coeffs)
529  type(sll_t_spline_pp_3d), intent(inout):: self
530  sll_int32, intent(in) :: n_cells(3)
531  sll_real64, intent(in) :: b_coeffs(n_cells(1)*n_cells(2), n_cells(2))
532  sll_real64, intent(out):: pp_coeffs((self%spline1%degree + 1)*(self%spline2%degree + 1)*(self%spline3%degree + 1), n_cells(1)*n_cells(2)*n_cells(3))
533  sll_int32 :: i, j, k
534 
535  do k = 1, n_cells(3)
536  do j = 1, n_cells(2)
537  do i = 1, n_cells(1)
538  call sll_s_spline_pp_b_to_pp_3d_cell(self, n_cells, b_coeffs, pp_coeffs, i, j, k)
539  end do
540  end do
541  end do
542 
543  end subroutine sll_s_spline_pp_b_to_pp_3d
544 
546  subroutine sll_s_spline_pp_b_to_pp_3d_cell(self, n_cells, b_coeffs, pp_coeffs, i, j, k)
547  type(sll_t_spline_pp_3d), intent(inout):: self
548  sll_int32, intent(in) :: n_cells(3)
549  sll_real64, intent(in) :: b_coeffs(n_cells(1)*n_cells(2)*n_cells(3))
550  sll_real64, intent(inout) :: pp_coeffs((self%spline1%degree + 1)*(self%spline2%degree + 1)*(self%spline3%degree + 1), n_cells(1)*n_cells(2)*n_cells(3))
551  sll_int32, intent(in) :: i, j, k
552  !local variables
553  sll_int32 :: l, m, n
554  sll_int32 :: degree1, degree2, degree3, degp1, degp2, degp3
555  degree1 = self%spline1%degree
556  degree2 = self%spline2%degree
557  degree3 = self%spline3%degree
558  degp1 = degree1 + 1
559  degp2 = degree2 + 1
560  degp3 = degree3 + 1
562  if (i > degree1) then
563  if (j > degree2) then
564  if (k > degree3) then
565  do m = 0, degree2
566  do n = 0, degree3
567  self%spline1%scratch_b = b_coeffs(i - degree1 + (j - degp2 + m)*n_cells(1) + (k - degp3 + n)*n_cells(1)*n_cells(2):i + (j - degp2 + m)*n_cells(1) + (k - degp3 + n)*n_cells(1)*n_cells(2))
568  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline1, self%spline1%scratch_b, self%scratch_pp(1 + m*degp1 + n*degp1*degp2:degp1*(1 + m + n*degp2)))
569  end do
570  end do
571  else
573  do m = 0, degree2
574  do n = 0, degree3
575  self%spline1%scratch_b = b_coeffs(i - degree1 + (j - degp2 + m)*n_cells(1) + modulo(k - degp3 + n, n_cells(3))*n_cells(1)*n_cells(2):i + (j - degp2 + m)*n_cells(1) + modulo(k - degp3 + n, n_cells(3))*n_cells(1)*n_cells(2))
576  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline1, self%spline1%scratch_b, self%scratch_pp(1 + m*degp1 + n*degp1*degp2:degp1*(1 + m + n*degp2)))
577  end do
578  end do
579  end if
580  else
582  do m = 0, degree2
583  do n = 0, degree3
584  self%spline1%scratch_b = b_coeffs(i - degree1 + modulo(j - degp2 + m, n_cells(2))*n_cells(1) + modulo(k - degp3 + n, n_cells(3))*n_cells(1)*n_cells(2):i + modulo(j - degp2 + m, n_cells(2))*n_cells(1) + modulo(k - degp3 + n, n_cells(3))*n_cells(1)*n_cells(2))
585  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline1, self%spline1%scratch_b, self%scratch_pp(1 + m*degp1 + n*degp1*degp2:degp1*(1 + m + n*degp2)))
586  end do
587  end do
588  end if
589  else
591  do m = 0, degree2
592  do n = 0, degree3
593  do l = 0, degree1
594  self%spline1%scratch_b(l + 1) = b_coeffs(modulo(i - degp1 + l, n_cells(1)) + 1 + modulo(j - degp2 + m, n_cells(2))*n_cells(1) + modulo(k - degp3 + n, n_cells(3))*n_cells(1)*n_cells(2))
595  end do
596  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline1, self%spline1%scratch_b, self%scratch_pp(1 + m*degp1 + n*degp1*degp2:degp1*(1 + m + n*degp2)))
597  end do
598  end do
599  end if
600 
602  do l = 1, degp1
603  do n = 1, degp3
604  do m = 1, degp2
605  self%spline2%scratch_b(m) = self%scratch_pp(l + degp1*(m - 1) + degp1*degp2*(n - 1))
606  end do
607  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline2, self%spline2%scratch_b, self%spline2%scratch_p)
608  do m = 1, degp2
609  self%scratch_pp(l + degp1*(m - 1) + degp1*degp2*(n - 1)) = self%spline2%scratch_p(m)
610  end do
611  end do
612  end do
613 
615  do l = 1, degp1
616  do m = 1, degp2
617  do n = 1, degp3
618  self%spline3%scratch_b(n) = self%scratch_pp(l + degp1*(m - 1) + degp1*degp2*(n - 1))
619  end do
620  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline3, self%spline3%scratch_b, self%spline3%scratch_p)
621  do n = 1, degp3
622  pp_coeffs(l + degp1*(m - 1) + degp1*degp2*(n - 1), i + n_cells(1)*(j - 1) + n_cells(1)*n_cells(2)*(k - 1)) = self%spline3%scratch_p(n)
623  end do
624  end do
625  end do
626 
627  end subroutine sll_s_spline_pp_b_to_pp_3d_cell
628 
630  subroutine sll_s_spline_pp_b_to_pp_3d_cella(self, n_cells, b_coeffs, pp_coeffs, i, j, k)
631  type(sll_t_spline_pp_3d), intent(inout):: self
632  sll_int32, intent(in) :: n_cells(3)
633  sll_real64, intent(in) :: b_coeffs((n_cells(1) + self%spline1%degree)*n_cells(2)*n_cells(3))
634  sll_real64, intent(inout) :: pp_coeffs((self%spline1%degree + 1)*(self%spline2%degree + 1)*(self%spline3%degree + 1), n_cells(1)*n_cells(2)*n_cells(3))
635  sll_int32, intent(in) :: i, j, k
636  !local variables
637  sll_int32 :: l, m, n
638  sll_int32 :: degree1, degree2, degree3, degp1, degp2, degp3, n_dofs
639  degree1 = self%spline1%degree
640  degree2 = self%spline2%degree
641  degree3 = self%spline3%degree
642  degp1 = degree1 + 1
643  degp2 = degree2 + 1
644  degp3 = degree3 + 1
645  n_dofs = n_cells(1) + degree1
647  if (i >= degree1 .and. i < n_cells(1) - degree1 + 2) then
648  if (j > degree2) then
649  if (k > degree3) then
650  do m = 0, degree2
651  do n = 0, degree3
652  self%spline1%scratch_b = b_coeffs(i + (j - degp2 + m)*n_dofs + (k - degp3 + n)*n_dofs*n_cells(2):i + degree1 + (j - degp2 + m)*n_dofs + (k - degp3 + n)*n_dofs*n_cells(2))
653  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline1, self%spline1%scratch_b, self%scratch_pp(1 + m*degp1 + n*degp1*degp2:degp1*(1 + m + n*degp2)))
654  end do
655  end do
656  else
658  do m = 0, degree2
659  do n = 0, degree3
660  self%spline1%scratch_b = b_coeffs(i + (j - degp2 + m)*n_dofs + modulo(k - degp3 + n, n_cells(3))*n_dofs*n_cells(2):i + degree1 + (j - degp2 + m)*n_dofs + modulo(k - degp3 + n, n_cells(3))*n_dofs*n_cells(2))
661  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline1, self%spline1%scratch_b, self%scratch_pp(1 + m*degp1 + n*degp1*degp2:degp1*(1 + m + n*degp2)))
662  end do
663  end do
664  end if
665  else
667  do m = 0, degree2
668  do n = 0, degree3
669  self%spline1%scratch_b = b_coeffs(i + modulo(j - degp2 + m, n_cells(2))*n_dofs + modulo(k - degp3 + n, n_cells(3))*n_dofs*n_cells(2):i + degree1 + modulo(j - degp2 + m, n_cells(2))*n_dofs + modulo(k - degp3 + n, n_cells(3))*n_dofs*n_cells(2))
670  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline1, self%spline1%scratch_b, self%scratch_pp(1 + m*degp1 + n*degp1*degp2:degp1*(1 + m + n*degp2)))
671  end do
672  end do
673  end if
674  else if (i < degree1) then
676  do m = 0, degree2
677  do n = 0, degree3
678  do l = 0, degree1
679  self%spline1%scratch_b(l + 1) = b_coeffs(i + l + modulo(j - degp2 + m, n_cells(2))*n_dofs + modulo(k - degp3 + n, n_cells(3))*n_dofs*n_cells(2))
680  end do
681  call sll_s_spline_pp_b_to_pp_1d_cella(degree1, self%spline1%poly_coeffs_boundary_left(:, :, i), self%spline1%scratch_b, self%scratch_pp(1 + m*degp1 + n*degp1*degp2:degp1*(1 + m + n*degp2)))
682  end do
683  end do
684  else if (i >= n_cells(1) - degree1 + 2) then
686  do m = 0, degree2
687  do n = 0, degree3
688  do l = 0, degree1
689  self%spline1%scratch_b(l + 1) = b_coeffs(i + l + modulo(j - degp2 + m, n_cells(2))*n_dofs + modulo(k - degp3 + n, n_cells(3))*n_dofs*n_cells(2))
690  end do
691  call sll_s_spline_pp_b_to_pp_1d_cella(degree1, self%spline1%poly_coeffs_boundary_right(:, :, i - (n_cells(1) - degree1 + 1)), self%spline1%scratch_b, self%scratch_pp(1 + m*degp1 + n*degp1*degp2:degp1*(1 + m + n*degp2)))
692  end do
693  end do
694  end if
695 
697  do l = 1, degp1
698  do n = 1, degp3
699  do m = 1, degp2
700  self%spline2%scratch_b(m) = self%scratch_pp(l + degp1*(m - 1) + degp1*degp2*(n - 1))
701  end do
702  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline2, self%spline2%scratch_b, self%spline2%scratch_p)
703  do m = 1, degp2
704  self%scratch_pp(l + degp1*(m - 1) + degp1*degp2*(n - 1)) = self%spline2%scratch_p(m)
705  end do
706  end do
707  end do
708 
710  do l = 1, degp1
711  do m = 1, degp2
712  do n = 1, degp3
713  self%spline3%scratch_b(n) = self%scratch_pp(l + degp1*(m - 1) + degp1*degp2*(n - 1))
714  end do
715  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline3, self%spline3%scratch_b, self%spline3%scratch_p)
716  do n = 1, degp3
717  pp_coeffs(l + degp1*(m - 1) + degp1*degp2*(n - 1), i + n_cells(1)*(j - 1) + n_cells(1)*n_cells(2)*(k - 1)) = self%spline3%scratch_p(n)
718  end do
719  end do
720  end do
721  end subroutine sll_s_spline_pp_b_to_pp_3d_cella
722 
724  subroutine sll_s_spline_pp_b_to_pp_3d_cella2f(self, n_cells, b_coeffs, pp_coeffs, i, j, k)
725  type(sll_t_spline_pp_3d), intent(inout):: self
726  sll_int32, intent(in) :: n_cells(3)
727  sll_real64, intent(in) :: b_coeffs((n_cells(1) + self%spline1%degree)*n_cells(2)*(n_cells(3) + self%spline3%degree))
728  sll_real64, intent(inout) :: pp_coeffs((self%spline1%degree + 1)*(self%spline2%degree + 1)*(self%spline3%degree + 1), n_cells(1)*n_cells(2)*n_cells(3))
729  sll_int32, intent(in) :: i, j, k
730  !local variables
731  sll_int32 :: l, m, n
732  sll_int32 :: deg(3), degp(3), n_dofs(3), upper(3)
733  deg(1) = self%spline1%degree
734  deg(2) = self%spline2%degree
735  deg(3) = self%spline3%degree
736  degp = deg + 1
737  n_dofs = n_cells + deg
738  n_dofs(2) = n_cells(2)
739  upper = n_cells - deg + 1
741  if (i >= deg(1) .and. i <= upper(1)) then
742  if (j > deg(2)) then
743  do n = 0, deg(3)
744  do m = 0, deg(2)
745  self%spline1%scratch_b = b_coeffs(i + (j - degp(2) + m)*n_dofs(1) + (k - 1 + n)*n_dofs(1)*n_dofs(2):i + deg(1) + (j - degp(2) + m)*n_dofs(1) + (k - 1 + n)*n_dofs(1)*n_dofs(2))
746  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline1, self%spline1%scratch_b, self%scratch_pp(1 + m*degp(1) + n*degp(1)*degp(2):degp(1)*(1 + m + n*degp(2))))
747  end do
748  end do
749  else
751  do n = 0, deg(3)
752  do m = 0, deg(2)
753  self%spline1%scratch_b = b_coeffs(i + modulo(j - degp(2) + m, n_cells(2))*n_dofs(1) + (k - 1 + n)*n_dofs(1)*n_dofs(2):i + deg(1) + modulo(j - degp(2) + m, n_cells(2))*n_dofs(1) + (k - 1 + n)*n_dofs(1)*n_dofs(2))
754  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline1, self%spline1%scratch_b, self%scratch_pp(1 + m*degp(1) + n*degp(1)*degp(2):degp(1)*(1 + m + n*degp(2))))
755  end do
756  end do
757  end if
758  else if (i < deg(1)) then
759  do n = 0, deg(3)
760  do m = 0, deg(2)
761  do l = 0, deg(1)
762  self%spline1%scratch_b(l + 1) = b_coeffs(i + l + modulo(j - degp(2) + m, n_cells(2))*n_dofs(1) + (k + n - 1)*n_dofs(1)*n_dofs(2))
763  end do
764  call sll_s_spline_pp_b_to_pp_1d_cella(deg(1), self%spline1%poly_coeffs_boundary_left(:, :, i), self%spline1%scratch_b, self%scratch_pp(1 + m*degp(1) + n*degp(1)*degp(2):degp(1)*(1 + m + n*degp(2))))
765  end do
766  end do
767  else if (i > upper(1)) then
768  do n = 0, deg(3)
769  do m = 0, deg(2)
770  do l = 0, deg(1)
771  self%spline1%scratch_b(l + 1) = b_coeffs(i + l + modulo(j - degp(2) + m, n_cells(2))*n_dofs(1) + (k + n - 1)*n_dofs(1)*n_dofs(2))
772  end do
773  call sll_s_spline_pp_b_to_pp_1d_cella(deg(1), self%spline1%poly_coeffs_boundary_right(:, :, i - upper(1)), self%spline1%scratch_b, self%scratch_pp(1 + m*degp(1) + n*degp(1)*degp(2):degp(1)*(1 + m + n*degp(2))))
774  end do
775  end do
776  end if
777 
779  do n = 0, deg(3)
780  do l = 0, deg(1)
781  do m = 0, deg(2)
782  self%spline2%scratch_b(m + 1) = self%scratch_pp(1 + l + degp(1)*m + degp(1)*degp(2)*n)
783  end do
784  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline2, self%spline2%scratch_b, self%spline2%scratch_p)
785  do m = 0, deg(2)
786  self%scratch_pp(1 + l + m*degp(1) + n*degp(1)*degp(2)) = self%spline2%scratch_p(m + 1)
787  end do
788  end do
789  end do
790 
792  if (k >= deg(3) .and. k <= upper(3)) then
793  do m = 0, deg(2)
794  do l = 0, deg(1)
795  do n = 0, deg(3)
796  self%spline3%scratch_b(n + 1) = self%scratch_pp(1 + l + degp(1)*m + degp(1)*degp(2)*n)
797  end do
798  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline3, self%spline3%scratch_b, self%spline3%scratch_p)
799  do n = 0, deg(3)
800  pp_coeffs(1 + l + m*degp(1) + n*degp(1)*degp(2), i + (j - 1)*n_cells(1) + (k - 1)*n_cells(1)*n_cells(2)) = self%spline3%scratch_p(n + 1)
801  end do
802  end do
803  end do
804  else if (k < deg(3)) then
805  do m = 0, deg(2)
806  do l = 0, deg(1)
807  do n = 0, deg(3)
808  self%spline3%scratch_b(n + 1) = self%scratch_pp(1 + l + degp(1)*m + degp(1)*degp(2)*n)
809  end do
810  call sll_s_spline_pp_b_to_pp_1d_cella(deg(3), self%spline3%poly_coeffs_boundary_left(:, :, k), self%spline3%scratch_b, self%spline3%scratch_p)
811  do n = 0, deg(3)
812  pp_coeffs(1 + l + m*degp(1) + n*degp(1)*degp(2), i + (j - 1)*n_cells(1) + (k - 1)*n_cells(1)*n_cells(2)) = self%spline3%scratch_p(n + 1)
813  end do
814  end do
815  end do
816  else if (k > upper(3)) then
817  do m = 0, deg(2)
818  do l = 0, deg(1)
819  do n = 0, deg(3)
820  self%spline3%scratch_b(n + 1) = self%scratch_pp(1 + l + degp(1)*m + degp(1)*degp(2)*n)
821  end do
822  call sll_s_spline_pp_b_to_pp_1d_cella(deg(3), self%spline3%poly_coeffs_boundary_right(:, :, k - upper(3)), self%spline3%scratch_b, self%spline3%scratch_p)
823  do n = 0, deg(3)
824  pp_coeffs(1 + l + m*degp(1) + n*degp(1)*degp(2), i + (j - 1)*n_cells(1) + (k - 1)*n_cells(1)*n_cells(2)) = self%spline3%scratch_p(n + 1)
825  end do
826  end do
827  end do
828  end if
829 
830  end subroutine sll_s_spline_pp_b_to_pp_3d_cella2f
831 
833  subroutine sll_s_spline_pp_b_to_pp_3d_cellaf(self, n_cells, b_coeffs, pp_coeffs, i, j, k)
834  type(sll_t_spline_pp_3d), intent(inout):: self
835  sll_int32, intent(in) :: n_cells(3)
836  sll_real64, intent(in) :: b_coeffs((n_cells(1) + self%spline1%degree)*(n_cells(2) + self%spline2%degree)*(n_cells(3) + self%spline3%degree))
837  sll_real64, intent(inout) :: pp_coeffs((self%spline1%degree + 1)*(self%spline2%degree + 1)*(self%spline3%degree + 1), n_cells(1)*n_cells(2)*n_cells(3))
838  sll_int32, intent(in) :: i, j, k
839  !local variables
840  sll_int32 :: l, m, n
841  sll_int32 :: deg(3), degp(3), n_dofs(3), upper(3)
842  deg(1) = self%spline1%degree
843  deg(2) = self%spline2%degree
844  deg(3) = self%spline3%degree
845  degp = deg + 1
846  n_dofs = n_cells + deg
847  upper = n_cells - deg + 1
849  if (i >= deg(1) .and. i <= upper(1)) then
850  do n = 0, deg(3)
851  do m = 0, deg(2)
852  self%spline1%scratch_b = b_coeffs(i + (j - 1 + m)*n_dofs(1) + (k - 1 + n)*n_dofs(1)*n_dofs(2):i + deg(1) + (j - 1 + m)*n_dofs(1) + (k - 1 + n)*n_dofs(1)*n_dofs(2))
853  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline1, self%spline1%scratch_b, self%scratch_pp(1 + m*degp(1) + n*degp(1)*degp(2):degp(1)*(1 + m + n*degp(2))))
854  end do
855  end do
856  else if (i < deg(1)) then
857  do n = 0, deg(3)
858  do m = 0, deg(2)
859  do l = 0, deg(1)
860  self%spline1%scratch_b(l + 1) = b_coeffs(i + l + (j + m - 1)*n_dofs(1) + (k + n - 1)*n_dofs(1)*n_dofs(2))
861  end do
862  call sll_s_spline_pp_b_to_pp_1d_cella(deg(1), self%spline1%poly_coeffs_boundary_left(:, :, i), self%spline1%scratch_b, self%scratch_pp(1 + m*degp(1) + n*degp(1)*degp(2):degp(1)*(1 + m + n*degp(2))))
863  end do
864  end do
865  else if (i > upper(1)) then
866  do n = 0, deg(3)
867  do m = 0, deg(2)
868  do l = 0, deg(1)
869  self%spline1%scratch_b(l + 1) = b_coeffs(i + l + (j + m - 1)*n_dofs(1) + (k + n - 1)*n_dofs(1)*n_dofs(2))
870  end do
871  call sll_s_spline_pp_b_to_pp_1d_cella(deg(1), self%spline1%poly_coeffs_boundary_right(:, :, i - upper(1)), self%spline1%scratch_b, self%scratch_pp(1 + m*degp(1) + n*degp(1)*degp(2):degp(1)*(1 + m + n*degp(2))))
872  end do
873  end do
874  end if
875 
877  if (j >= deg(2) .and. j <= upper(2)) then
878  do n = 0, deg(3)
879  do l = 0, deg(1)
880  do m = 0, deg(2)
881  self%spline2%scratch_b(m + 1) = self%scratch_pp(1 + l + degp(1)*m + degp(1)*degp(2)*n)
882  end do
883  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline2, self%spline2%scratch_b, self%spline2%scratch_p)
884  do m = 0, deg(2)
885  self%scratch_pp(1 + l + m*degp(1) + n*degp(1)*degp(2)) = self%spline2%scratch_p(m + 1)
886  end do
887  end do
888  end do
889  else if (j < deg(2)) then
890  do n = 0, deg(3)
891  do l = 0, deg(1)
892  do m = 0, deg(2)
893  self%spline2%scratch_b(m + 1) = self%scratch_pp(1 + l + degp(1)*m + degp(1)*degp(2)*n)
894  end do
895  call sll_s_spline_pp_b_to_pp_1d_cella(deg(2), self%spline2%poly_coeffs_boundary_left(:, :, j), self%spline2%scratch_b, self%spline2%scratch_p)
896  do m = 0, deg(2)
897  self%scratch_pp(1 + l + m*degp(1) + n*degp(1)*degp(2)) = self%spline2%scratch_p(m + 1)
898  end do
899  end do
900  end do
901  else if (j > upper(2)) then
902  do n = 0, deg(3)
903  do l = 0, deg(1)
904  do m = 0, deg(2)
905  self%spline2%scratch_b(m + 1) = self%scratch_pp(1 + l + degp(1)*m + degp(1)*degp(2)*n)
906  end do
907  call sll_s_spline_pp_b_to_pp_1d_cella(deg(2), self%spline2%poly_coeffs_boundary_right(:, :, j - upper(2)), self%spline2%scratch_b, self%spline2%scratch_p)
908  do m = 0, deg(2)
909  self%scratch_pp(1 + l + m*degp(1) + n*degp(1)*degp(2)) = self%spline2%scratch_p(m + 1)
910  end do
911  end do
912  end do
913  end if
914 
916  if (k >= deg(3) .and. k <= upper(3)) then
917  do m = 0, deg(2)
918  do l = 0, deg(1)
919  do n = 0, deg(3)
920  self%spline3%scratch_b(n + 1) = self%scratch_pp(1 + l + degp(1)*m + degp(1)*degp(2)*n)
921  end do
922  call sll_s_spline_pp_b_to_pp_1d_cell(self%spline3, self%spline3%scratch_b, self%spline3%scratch_p)
923  do n = 0, deg(3)
924  pp_coeffs(1 + l + m*degp(1) + n*degp(1)*degp(2), i + (j - 1)*n_cells(1) + (k - 1)*n_cells(1)*n_cells(2)) = self%spline3%scratch_p(n + 1)
925  end do
926  end do
927  end do
928  else if (k < deg(3)) then
929  do m = 0, deg(2)
930  do l = 0, deg(1)
931  do n = 0, deg(3)
932  self%spline3%scratch_b(n + 1) = self%scratch_pp(1 + l + degp(1)*m + degp(1)*degp(2)*n)
933  end do
934  call sll_s_spline_pp_b_to_pp_1d_cella(deg(3), self%spline3%poly_coeffs_boundary_left(:, :, k), self%spline3%scratch_b, self%spline3%scratch_p)
935  do n = 0, deg(3)
936  pp_coeffs(1 + l + m*degp(1) + n*degp(1)*degp(2), i + (j - 1)*n_cells(1) + (k - 1)*n_cells(1)*n_cells(2)) = self%spline3%scratch_p(n + 1)
937  end do
938  end do
939  end do
940  else if (k > upper(3)) then
941  do m = 0, deg(2)
942  do l = 0, deg(1)
943  do n = 0, deg(3)
944  self%spline3%scratch_b(n + 1) = self%scratch_pp(1 + l + degp(1)*m + degp(1)*degp(2)*n)
945  end do
946  call sll_s_spline_pp_b_to_pp_1d_cella(deg(3), self%spline3%poly_coeffs_boundary_right(:, :, k - upper(3)), self%spline3%scratch_b, self%spline3%scratch_p)
947  do n = 0, deg(3)
948  pp_coeffs(1 + l + m*degp(1) + n*degp(1)*degp(2), i + (j - 1)*n_cells(1) + (k - 1)*n_cells(1)*n_cells(2)) = self%spline3%scratch_p(n + 1)
949  end do
950  end do
951  end do
952  end if
953 
954  end subroutine sll_s_spline_pp_b_to_pp_3d_cellaf
955 
957  subroutine sll_s_spline_pp_b_to_pp_3d_clamped(self, n_cells, b_coeffs, pp_coeffs)
958  type(sll_t_spline_pp_3d), intent(inout):: self
959  sll_int32, intent(in) :: n_cells(3)
960  sll_real64, intent(in) :: b_coeffs((n_cells(1) + self%spline1%degree)*n_cells(2)*n_cells(3))
961  sll_real64, intent(out):: pp_coeffs((self%spline1%degree + 1)*(self%spline2%degree + 1)*(self%spline3%degree + 1), n_cells(1)*n_cells(2)*n_cells(3))
962  sll_int32 :: i, j, k
963 
964  do k = 1, n_cells(3)
965  do j = 1, n_cells(2)
966  do i = 1, n_cells(1)
967  call sll_s_spline_pp_b_to_pp_3d_cella(self, n_cells, b_coeffs, pp_coeffs, i, j, k)
968  end do
969  end do
970  end do
971  end subroutine sll_s_spline_pp_b_to_pp_3d_clamped
972 
974  subroutine sll_s_spline_pp_b_to_pp_3d_clamped_2full(self, n_cells, b_coeffs, pp_coeffs)
975  type(sll_t_spline_pp_3d), intent(inout):: self
976  sll_int32, intent(in) :: n_cells(3)
977  sll_real64, intent(in) :: b_coeffs((n_cells(1) + self%spline1%degree)*n_cells(2)*(n_cells(3) + self%spline3%degree))
978  sll_real64, intent(out):: pp_coeffs((self%spline1%degree + 1)*(self%spline2%degree + 1)*(self%spline3%degree + 1), n_cells(1)*n_cells(2)*n_cells(3))
979  sll_int32 :: i, j, k
980 
981  do k = 1, n_cells(3)
982  do j = 1, n_cells(2)
983  do i = 1, n_cells(1)
984  call sll_s_spline_pp_b_to_pp_3d_cella2f(self, n_cells, b_coeffs, pp_coeffs, i, j, k)
985  end do
986  end do
987  end do
988 
989  end subroutine sll_s_spline_pp_b_to_pp_3d_clamped_2full
990 
992  subroutine sll_s_spline_pp_b_to_pp_3d_clamped_full(self, n_cells, b_coeffs, pp_coeffs)
993  type(sll_t_spline_pp_3d), intent(inout):: self
994  sll_int32, intent(in) :: n_cells(3)
995  sll_real64, intent(in) :: b_coeffs((n_cells(1) + self%spline1%degree)*(n_cells(2) + self%spline2%degree)*(n_cells(3) + self%spline3%degree))
996  sll_real64, intent(out):: pp_coeffs((self%spline1%degree + 1)*(self%spline2%degree + 1)*(self%spline3%degree + 1), n_cells(1)*n_cells(2)*n_cells(3))
997  sll_int32 :: i, j, k
998 
999  do k = 1, n_cells(3)
1000  do j = 1, n_cells(2)
1001  do i = 1, n_cells(1)
1002  call sll_s_spline_pp_b_to_pp_3d_cellaf(self, n_cells, b_coeffs, pp_coeffs, i, j, k)
1003  end do
1004  end do
1005  end do
1006 
1007  end subroutine sll_s_spline_pp_b_to_pp_3d_clamped_full
1008 
1009  subroutine sll_s_spline_evaluate_basis_b_form_1d_clamped(self, n_cells, b_coeffs, val)
1010  type(sll_t_spline_pp_1d), intent(in):: self
1011  sll_int32, intent(in) :: n_cells
1012  sll_real64, intent(in) :: b_coeffs(n_cells + self%degree)
1013  sll_real64, intent(out):: val(:)
1014  !local variables
1015  sll_real64 :: pp_coeffs(self%degree + 1, n_cells)
1016 
1017  call sll_s_spline_pp_b_to_pp_1d_clamped(self, n_cells, b_coeffs, pp_coeffs)
1018 
1019  call sll_s_spline_evaluate_basis_pp_form_1d(self, n_cells, pp_coeffs, val)
1020 
1021  end subroutine sll_s_spline_evaluate_basis_b_form_1d_clamped
1022 
1023  subroutine sll_s_spline_evaluate_basis_b_form_1d_periodic(self, n_cells, b_coeffs, val)
1024  type(sll_t_spline_pp_1d), intent(in):: self
1025  sll_int32, intent(in) :: n_cells
1026  sll_real64, intent(in) :: b_coeffs(n_cells)
1027  sll_real64, intent(out):: val(:)
1028  !local variables
1029  sll_real64 :: pp_coeffs(self%degree + 1, n_cells)
1030 
1031  call sll_s_spline_pp_b_to_pp_1d(self, n_cells, b_coeffs, pp_coeffs)
1032 
1033  call sll_s_spline_evaluate_basis_pp_form_1d(self, n_cells, pp_coeffs, val)
1034 
1035  end subroutine sll_s_spline_evaluate_basis_b_form_1d_periodic
1036 
1037  subroutine sll_s_spline_evaluate_basis_pp_form_1d(self, n_cells, pp_coeffs, val)
1038  type(sll_t_spline_pp_1d), intent(in):: self
1039  sll_int32, intent(in) :: n_cells
1040  sll_real64, intent(in) :: pp_coeffs(self%degree + 1, n_cells)
1041  sll_real64, intent(out):: val(:)
1042  !local variables
1043  sll_int32 :: i
1044 
1045  do i = 1, n_cells
1046  val(i) = pp_coeffs(self%degree + 1, i)
1047  !val(i)=sll_f_spline_pp_horner_1d(self%degree, pp_coeffs, 0._f64, i)
1048  end do
1049  val(n_cells + 1) = sll_f_spline_pp_horner_1d(self%degree, pp_coeffs, 1._f64, n_cells)
1050 
1051  end subroutine sll_s_spline_evaluate_basis_pp_form_1d
1052 
1053  subroutine sll_s_spline_evaluate_basis_b_form_3d_clamped(self, n_cells, b_coeffs, val)
1054  type(sll_t_spline_pp_3d), intent(inout):: self
1055  sll_int32, intent(in) :: n_cells(3)
1056  sll_real64, intent(in) :: b_coeffs((n_cells(1) + self%spline1%degree)*n_cells(2)*n_cells(3))
1057  sll_real64, intent(out):: val(:)
1058 
1059  call sll_s_spline_pp_b_to_pp_3d_clamped(self, n_cells, b_coeffs, self%scratch_coeffs)
1060 
1061  call sll_s_spline_evaluate_basis_pp_form_3d(self, n_cells, self%scratch_coeffs, val)
1062 
1063  end subroutine sll_s_spline_evaluate_basis_b_form_3d_clamped
1064 
1065  subroutine sll_s_spline_evaluate_basis_b_form_3d_periodic(self, n_cells, b_coeffs, val)
1066  type(sll_t_spline_pp_3d), intent(inout):: self
1067  sll_int32, intent(in) :: n_cells(3)
1068  sll_real64, intent(in) :: b_coeffs(n_cells(1)*n_cells(2)*n_cells(3))
1069  sll_real64, intent(out):: val(:)
1070 
1071  call sll_s_spline_pp_b_to_pp_3d(self, n_cells, b_coeffs, self%scratch_coeffs)
1072 
1073  call sll_s_spline_evaluate_basis_pp_form_3d(self, n_cells, self%scratch_coeffs, val)
1074 
1075  end subroutine sll_s_spline_evaluate_basis_b_form_3d_periodic
1076 
1077  subroutine sll_s_spline_evaluate_basis_pp_form_3d(self, n_cells, pp_coeffs, val)
1078  type(sll_t_spline_pp_3d), intent(in):: self
1079  sll_int32, intent(in) :: n_cells(3)
1080  sll_real64, intent(in) :: pp_coeffs((self%spline1%degree + 1)*(self%spline2%degree + 1)*(self%spline3%degree + 1), n_cells(1)*n_cells(2)*n_cells(3))
1081  sll_real64, intent(out):: val(:)
1082  !local variables
1083  sll_int32 :: i, j, k
1084 
1085  do k = 1, n_cells(3)
1086  do j = 1, n_cells(2)
1087  do i = 1, n_cells(1)
1088  !val(i+(j-1)*n_cells(1)+(k-1)*n_cells(1)*n_cells(2)) = pp_coeffs(self%spline1%degree+1,i+(j-1)*n_cells(1)+(k-1)*n_cells(1)*n_cells(2))
1089  val(i + (j - 1)*(n_cells(1) + 1) + (k - 1)*(n_cells(1) + 1)*(n_cells(2) + 1)) = sll_f_spline_pp_horner_3d([self%spline1%degree, self%spline2%degree, self%spline3%degree], pp_coeffs, [0._f64, 0._f64, 0._f64], [i, j, k], n_cells)
1090  end do
1091  val(n_cells(1) + 1 + (j - 1)*(n_cells(1) + 1) + (k - 1)*(n_cells(1) + 1)*(n_cells(2) + 1)) = sll_f_spline_pp_horner_3d([self%spline1%degree, self%spline2%degree, self%spline3%degree], pp_coeffs, [1._f64, 0._f64, 0._f64], [n_cells(1), j, k], n_cells)
1092  end do
1093  val(n_cells(1) + 1 + n_cells(2)*(n_cells(1) + 1) + (k - 1)*(n_cells(1) + 1)*(n_cells(2) + 1)) = sll_f_spline_pp_horner_3d([self%spline1%degree, self%spline2%degree, self%spline3%degree], pp_coeffs, [1._f64, 1._f64, 0._f64], [n_cells(1), n_cells(2), k], n_cells)
1094  end do
1095  val(n_cells(1) + 1 + n_cells(2)*(n_cells(1) + 1) + n_cells(3)*(n_cells(1) + 1)*(n_cells(2) + 1)) = sll_f_spline_pp_horner_3d([self%spline1%degree, self%spline2%degree, self%spline3%degree], pp_coeffs, [1._f64, 1._f64, 1._f64], [n_cells(1), n_cells(2), n_cells(3)], n_cells)
1096 
1097  do i = 1, n_cells(1)
1098  do k = 1, n_cells(3)
1099  val(i + n_cells(2)*(n_cells(1) + 1) + (k - 1)*(n_cells(1) + 1)*(n_cells(2) + 1)) = sll_f_spline_pp_horner_3d([self%spline1%degree, self%spline2%degree, self%spline3%degree], pp_coeffs, [0._f64, 1._f64, 0._f64], [i, n_cells(2), k], n_cells)
1100  end do
1101  val(i + n_cells(2)*(n_cells(1) + 1) + n_cells(3)*(n_cells(1) + 1)*(n_cells(2) + 1)) = sll_f_spline_pp_horner_3d([self%spline1%degree, self%spline2%degree, self%spline3%degree], pp_coeffs, [0._f64, 1._f64, 1._f64], [i, n_cells(2), n_cells(3)], n_cells)
1102  end do
1103  do j = 1, n_cells(2)
1104  do i = 1, n_cells(1)
1105  val(i + (j - 1)*(n_cells(1) + 1) + n_cells(3)*(n_cells(1) + 1)*(n_cells(2) + 1)) = sll_f_spline_pp_horner_3d([self%spline1%degree, self%spline2%degree, self%spline3%degree], pp_coeffs, [0._f64, 0._f64, 1._f64], [i, j, n_cells(3)], n_cells)
1106  end do
1107  val(n_cells(1) + 1 + (j - 1)*(n_cells(1) + 1) + n_cells(3)*(n_cells(1) + 1)*(n_cells(2) + 1)) = sll_f_spline_pp_horner_3d([self%spline1%degree, self%spline2%degree, self%spline3%degree], pp_coeffs, [1._f64, 0._f64, 1._f64], [n_cells(1), j, n_cells(3)], n_cells)
1108  end do
1109 
1110  end subroutine sll_s_spline_evaluate_basis_pp_form_3d
1111 
1113  subroutine sll_s_spline_pp_horner_m_1d(self, val, degree, x)
1114  type(sll_t_spline_pp_1d), intent(in):: self
1115  sll_real64, intent(out):: val(:)
1116  sll_int32, intent(in) :: degree
1117  sll_real64, intent(in) :: x
1118  sll_int32 :: i
1119 
1120  do i = 1, degree + 1!size(val)
1121  val(i) = sll_f_spline_pp_horner_1d(degree, self%poly_coeffs, x, i)
1122  end do
1123  end subroutine sll_s_spline_pp_horner_m_1d
1124 
1126  subroutine sll_s_spline_pp_horner_m_2d(self, val, degree, x)
1127  type(sll_t_spline_pp_2d), intent(in):: self
1128  sll_real64, intent(out):: val(:, :)
1129  sll_int32, intent(in) :: degree(2)
1130  sll_real64, intent(in) :: x(2)
1131  sll_int32 :: i
1132 
1133  do i = 1, degree(1) + 1!size(val,1)
1134  val(i, 1) = sll_f_spline_pp_horner_1d(degree(1), self%spline1%poly_coeffs, x(1), i)
1135  end do
1136  do i = 1, degree(2) + 1!size(val,2)
1137  val(i, 2) = sll_f_spline_pp_horner_1d(degree(2), self%spline2%poly_coeffs, x(2), i)
1138  end do
1139  end subroutine sll_s_spline_pp_horner_m_2d
1140 
1142  subroutine sll_s_spline_pp_horner_m_3d(self, val, degree, x)
1143  type(sll_t_spline_pp_3d), intent(in):: self
1144  sll_real64, intent(out):: val(:, :)
1145  sll_int32, intent(in) :: degree(3)
1146  sll_real64, intent(in) :: x(3)
1147  sll_int32 :: i
1148 
1149  do i = 1, degree(1) + 1!size(val,1)
1150  val(i, 1) = sll_f_spline_pp_horner_1d(degree(1), self%spline1%poly_coeffs, x(1), i)
1151  end do
1152  do i = 1, degree(2) + 1!size(val,2)
1153  val(i, 2) = sll_f_spline_pp_horner_1d(degree(2), self%spline2%poly_coeffs, x(2), i)
1154  end do
1155  do i = 1, degree(3) + 1!size(val,3)
1156  val(i, 3) = sll_f_spline_pp_horner_1d(degree(3), self%spline3%poly_coeffs, x(3), i)
1157  end do
1158  end subroutine sll_s_spline_pp_horner_m_3d
1159 
1161  subroutine sll_s_spline_pp_horner_primitive_1d(val, degree, pp_coeffs, x)
1162  sll_real64, intent(out):: val(:)
1163  sll_int32, intent(in) :: degree
1164  sll_real64, intent(in) :: pp_coeffs(:, :)
1165  sll_real64, intent(in) :: x
1166  sll_int32 :: i
1167 
1168  do i = 1, size(val)
1169  val(i) = sll_f_spline_pp_horner_1d(degree, pp_coeffs, x, i)*x
1170  end do
1171  end subroutine sll_s_spline_pp_horner_primitive_1d
1172 
1174  function sll_f_spline_pp_horner_1d(degree, pp_coeffs, x, index) result(res)
1175  sll_int32, intent(in) :: degree
1176  sll_real64, intent(in) :: pp_coeffs(:, :)
1177  sll_real64, intent(in) :: x
1178  sll_int32, intent(in) :: index
1179  sll_real64 :: res
1180  sll_int32 :: i
1181 
1182  !print*, pp_coeffs(:,index)
1183  !print*, x
1184  res = pp_coeffs(1, index)
1185  do i = 1, degree
1186  res = res*x + pp_coeffs(i + 1, index)
1187  end do
1188  end function sll_f_spline_pp_horner_1d
1189 
1191  function sll_f_spline_pp_horner_derivative_1d(degree, pp_coeffs, x, index) result(res)
1192  sll_int32, intent(in) :: degree
1193  sll_real64, intent(in) :: pp_coeffs(:, :)
1194  sll_real64, intent(in) :: x
1195  sll_int32, intent(in) :: index
1196  sll_real64 :: res
1197  sll_int32 :: i
1198 
1199  !print*, pp_coeffs(:,index)
1200  !print*, x
1201  res = pp_coeffs(1, index)*real(degree, f64)
1202  do i = 1, degree - 1
1203  res = res*x + pp_coeffs(i + 1, index)*real(degree - i, f64)
1204  end do
1205  end function sll_f_spline_pp_horner_derivative_1d
1206 
1208  function sll_f_spline_pp_horner_2d(degree, pp_coeffs, x, indices, n_cells) result(res)
1209  sll_int32, intent(in) :: degree(2)
1210  sll_real64, intent(in) :: pp_coeffs(:, :)
1211  sll_real64, intent(in) :: x(2)
1212  sll_int32, intent(in) :: indices(2)
1213  sll_int32, intent(in) :: n_cells(2)
1214  sll_real64 :: res
1215  sll_real64 :: pp_coeffs_1d(degree(2) + 1, 1)
1216  sll_int32 :: i
1218  do i = 0, degree(2)
1219  pp_coeffs_1d(i + 1, 1) = sll_f_spline_pp_horner_1d(degree(1), pp_coeffs(1 + i*(degree(1) + 1):(degree(1) + 1)*(i + 1), 1 + (indices(2) - 1)*n_cells(1):n_cells(1)*indices(2)), x(1), indices(1))
1220  end do
1222  res = sll_f_spline_pp_horner_1d(degree(2), pp_coeffs_1d, x(2), 1)
1223  end function sll_f_spline_pp_horner_2d
1224 
1226  function sll_f_spline_pp_horner_3d(degree, pp_coeffs, x, indices, n_cells) result(res)
1227  sll_int32, intent(in) :: degree(3)
1228  sll_real64, intent(in) :: pp_coeffs(:, :)
1229  sll_real64, intent(in) :: x(3)
1230  sll_int32, intent(in) :: indices(3)
1231  sll_int32, intent(in) :: n_cells(3)
1232  sll_real64 :: res
1233  sll_real64 :: pp_coeffs_2d((degree(2) + 1)*(degree(3) + 1), 1)
1234  sll_int32 :: i, j
1235  sll_int32 :: degp1, degp2
1236  degp1 = degree(1) + 1
1237  degp2 = degree(2) + 1
1239  do j = 0, degree(3)
1240  do i = 0, degree(2)
1241  pp_coeffs_2d(1 + i + j*degp2, 1) = sll_f_spline_pp_horner_1d(degree(1), pp_coeffs(1 + i*degp1 + j*degp1*degp2:degp1 + i*degp1 + j*degp1*degp2, 1 + n_cells(1)*(indices(2) - 1 + (indices(3) - 1)*n_cells(2)):n_cells(1)*(indices(2) + (indices(3) - 1)*n_cells(2))), x(1), indices(1))
1242  end do
1243  end do
1245  res = sll_f_spline_pp_horner_2d(degree(2:3), pp_coeffs_2d, x(2:3), [1, 1], [1, 1])
1246 
1247  end function sll_f_spline_pp_horner_3d
1248 
1250  function sll_f_spline_pp_horner_3d_d1(degree, pp_coeffs, x, indices, n_cells) result(res)
1251  sll_int32, intent(in) :: degree(3)
1252  sll_real64, intent(in) :: pp_coeffs(:, :)
1253  sll_real64, intent(in) :: x(3)
1254  sll_int32, intent(in) :: indices(3)
1255  sll_int32, intent(in) :: n_cells(3)
1256  sll_real64 :: res
1257  sll_real64 :: pp_coeffs_2d((degree(2) + 1)*(degree(3) + 1), 1)
1258  sll_real64 :: pp_coeffs_1d(degree(3) + 1, 1)
1259  sll_int32 :: i, j
1260  sll_int32 :: degp1, degp2
1261  degp1 = degree(1) + 1
1262  degp2 = degree(2) + 1
1264  do j = 0, degree(3)
1265  do i = 0, degree(2)
1266  pp_coeffs_2d(1 + i + j*degp2, 1) = sll_f_spline_pp_horner_derivative_1d(degree(1), pp_coeffs(1 + i*degp1 + j*degp1*degp2:degp1 + i*degp1 + j*degp1*degp2, indices(1) + n_cells(1)*(indices(2) - 1 + (indices(3) - 1)*n_cells(2)):indices(1) + n_cells(1)*(indices(2) - 1 + (indices(3) - 1)*n_cells(2))), x(1), 1)
1267  end do
1268  end do
1270  res = sll_f_spline_pp_horner_2d(degree(2:3), pp_coeffs_2d, x(2:3), [1, 1], [1, 1])
1271 
1273  do i = 0, degree(3)
1274  pp_coeffs_1d(i + 1, 1) = sll_f_spline_pp_horner_1d(degree(2), pp_coeffs_2d(1 + i*(degree(2) + 1):(degree(2) + 1)*(i + 1), :), x(2), 1)
1275  end do
1277  res = sll_f_spline_pp_horner_1d(degree(3), pp_coeffs_1d, x(3), 1)*real(n_cells(1), f64)
1278 
1279  end function sll_f_spline_pp_horner_3d_d1
1280 
1282  function sll_f_spline_pp_horner_3d_d2(degree, pp_coeffs, x, indices, n_cells) result(res)
1283  sll_int32, intent(in) :: degree(3)
1284  sll_real64, intent(in) :: pp_coeffs(:, :)
1285  sll_real64, intent(in) :: x(3)
1286  sll_int32, intent(in) :: indices(3)
1287  sll_int32, intent(in) :: n_cells(3)
1288  sll_real64 :: res
1289  sll_real64 :: pp_coeffs_2d((degree(2) + 1)*(degree(3) + 1), 1)
1290  sll_real64 :: pp_coeffs_1d(degree(3) + 1, 1)
1291  sll_int32 :: i, j
1292  sll_int32 :: degp1, degp2
1293  degp1 = degree(1) + 1
1294  degp2 = degree(2) + 1
1296  do j = 0, degree(3)
1297  do i = 0, degree(2)
1298  pp_coeffs_2d(1 + i + j*degp2, 1) = sll_f_spline_pp_horner_1d(degree(1), pp_coeffs(1 + i*degp1 + j*degp1*degp2:degp1 + i*degp1 + j*degp1*degp2, indices(1) + n_cells(1)*(indices(2) - 1 + (indices(3) - 1)*n_cells(2)):indices(1) + n_cells(1)*(indices(2) - 1 + (indices(3) - 1)*n_cells(2))), x(1), 1)
1299  end do
1300  end do
1302  res = sll_f_spline_pp_horner_2d(degree(2:3), pp_coeffs_2d, x(2:3), [1, 1], [1, 1])
1303 
1305  do i = 0, degree(3)
1306  pp_coeffs_1d(i + 1, 1) = sll_f_spline_pp_horner_derivative_1d(degree(2), pp_coeffs_2d(1 + i*(degree(2) + 1):(degree(2) + 1)*(i + 1), :), x(2), 1)
1307  end do
1309  res = sll_f_spline_pp_horner_1d(degree(3), pp_coeffs_1d, x(3), 1)*real(n_cells(2), f64)
1310 
1311  end function sll_f_spline_pp_horner_3d_d2
1312 
1314  function sll_f_spline_pp_horner_3d_d3(degree, pp_coeffs, x, indices, n_cells) result(res)
1315  sll_int32, intent(in) :: degree(3)
1316  sll_real64, intent(in) :: pp_coeffs(:, :)
1317  sll_real64, intent(in) :: x(3)
1318  sll_int32, intent(in) :: indices(3)
1319  sll_int32, intent(in) :: n_cells(3)
1320  sll_real64 :: res
1321  sll_real64 :: pp_coeffs_2d((degree(2) + 1)*(degree(3) + 1), 1)
1322  sll_real64 :: pp_coeffs_1d(degree(3) + 1, 1)
1323  sll_int32 :: i, j
1324  sll_int32 :: degp1, degp2
1325  degp1 = degree(1) + 1
1326  degp2 = degree(2) + 1
1328  do j = 0, degree(3)
1329  do i = 0, degree(2)
1330  pp_coeffs_2d(1 + i + j*degp2, 1) = sll_f_spline_pp_horner_1d(degree(1), pp_coeffs(1 + i*degp1 + j*degp1*degp2:degp1 + i*degp1 + j*degp1*degp2, indices(1) + n_cells(1)*(indices(2) - 1 + (indices(3) - 1)*n_cells(2)):indices(1) + n_cells(1)*(indices(2) - 1 + (indices(3) - 1)*n_cells(2))), x(1), 1)
1331  end do
1332  end do
1334  res = sll_f_spline_pp_horner_2d(degree(2:3), pp_coeffs_2d, x(2:3), [1, 1], [1, 1])
1335 
1337  do i = 0, degree(3)
1338  pp_coeffs_1d(i + 1, 1) = sll_f_spline_pp_horner_1d(degree(2), pp_coeffs_2d(1 + i*(degree(2) + 1):(degree(2) + 1)*(i + 1), :), x(2), 1)
1339  end do
1341  res = sll_f_spline_pp_horner_derivative_1d(degree(3), pp_coeffs_1d, x(3), 1)*real(n_cells(3), f64)
1342 
1343  end function sll_f_spline_pp_horner_3d_d3
1344 
1346  function sll_f_spline_pp_horner_derivative_3d(degree, pp_coeffs, x, indices, n_cells, component) result(res)
1347  sll_int32, intent(in) :: degree(3)
1348  sll_real64, intent(in) :: pp_coeffs(:, :)
1349  sll_real64, intent(in) :: x(3)
1350  sll_int32, intent(in) :: indices(3)
1351  sll_int32, intent(in) :: n_cells(3)
1352  sll_int32, intent(in) :: component
1353  sll_real64 :: res
1354  sll_real64, allocatable :: coeffs(:, :), pp_coeffs_2d(:, :)
1355  sll_real64 :: exponents(maxval(degree) + 1, 3)
1356  sll_int32 :: h, i, j
1357  sll_int32 :: deg(3), degp1, degp2
1358  degp1 = degree(1) + 1
1359  degp2 = degree(2) + 1
1360 
1361  deg = degree
1362  deg(component) = degree(component) - 1
1363  allocate (coeffs((deg(1) + 1), 1))
1364  allocate (pp_coeffs_2d((deg(2) + 1)*(deg(3) + 1), 1))
1365 
1366  exponents = 1._f64
1367  do i = 1, degree(component)
1368  exponents(i, component) = real(degree(component) + 1 - i, f64)
1369  end do
1371  do j = 0, deg(3)
1372  do i = 0, deg(2)
1373  do h = 0, deg(1)
1374  coeffs(h + 1, 1) = pp_coeffs(1 + h + i*degp1 + j*degp1*degp2, indices(1) + (indices(2) - 1)*n_cells(1) + (indices(3) - 1)*n_cells(1)*n_cells(2))*exponents(j + 1, 3)*exponents(i + 1, 2)*exponents(h + 1, 1)
1375  end do
1376  pp_coeffs_2d(1 + i + j*(deg(2) + 1), 1) = sll_f_spline_pp_horner_1d(deg(1), coeffs, x(1), 1)
1377  end do
1378  end do
1379 
1381  res = sll_f_spline_pp_horner_2d(deg(2:3), pp_coeffs_2d, x(2:3), [1, 1], [1, 1])*real(n_cells(component), f64)
1382 
1383  end function sll_f_spline_pp_horner_derivative_3d
1384 
1386  subroutine sll_s_spline_pp_init_1d(self, degree, n_cells, boundary)
1387  type(sll_t_spline_pp_1d), intent(out) :: self
1388  sll_int32, intent(in) :: degree
1389  sll_int32, intent(in) :: n_cells
1390  sll_int32, intent(in), optional :: boundary
1391  sll_int32 :: ierr
1392 
1393  if (present(boundary)) then
1394  self%boundary_conditions = boundary
1395  end if
1396 
1397  sll_assert(n_cells >= degree)
1398  self%degree = degree
1399  self%n_cells = n_cells
1400  sll_allocate(self%poly_coeffs(degree + 1, degree + 1), ierr)
1401  sll_assert(ierr == 0)
1402  sll_allocate(self%poly_coeffs_fp(degree + 1, degree + 1), ierr)
1403  sll_assert(ierr == 0)
1404  sll_allocate(self%poly_coeffs_fpa(degree + 2, degree + 1), ierr)
1405  sll_assert(ierr == 0)
1406  sll_allocate(self%poly_coeffs_fd(degree, degree + 1), ierr)
1407  sll_assert(ierr == 0)
1408  sll_allocate(self%scratch_b(degree + 1), ierr)
1409  sll_assert(ierr == 0)
1410  sll_allocate(self%scratch_p(degree + 1), ierr)
1411  sll_assert(ierr == 0)
1412 
1413  sll_allocate(self%poly_coeffs_boundary_left(degree + 1, degree + 1, degree - 1), ierr)
1414  sll_assert(ierr == 0)
1415  sll_allocate(self%poly_coeffs_boundary_right(degree + 1, degree + 1, degree - 1), ierr)
1416  sll_assert(ierr == 0)
1417 
1418  sll_allocate(self%poly_coeffs_fd_boundary_left(degree, degree + 1, degree - 1), ierr)
1419  sll_assert(ierr == 0)
1420  sll_allocate(self%poly_coeffs_fd_boundary_right(degree, degree + 1, degree - 1), ierr)
1421  sll_assert(ierr == 0)
1422 
1423  select case (self%boundary_conditions)
1424  case (sll_p_boundary_periodic)
1425  self%n_coeffs = n_cells
1426  case (sll_p_boundary_clamped)
1427  self%n_coeffs = n_cells + degree
1428  case (sll_p_boundary_clamped_square)
1429  self%n_coeffs = n_cells + degree
1430  case (sll_p_boundary_clamped_cubic)
1431  self%n_coeffs = n_cells + degree
1432  case (sll_p_boundary_clamped_clampeddiri)
1433  self%n_coeffs = n_cells + degree - 1
1434  case (sll_p_boundary_clampeddiri)
1435  self%n_coeffs = n_cells + degree - 2
1436  case (sll_p_boundary_clampeddiri_clamped)
1437  self%n_coeffs = n_cells + degree - 1
1438  case default
1439  sll_error('sll_s_spline_pp_init_1d', 'Wrong boundary condition.')
1440  end select
1441 
1442  select case (self%degree)
1443  !poly_coefficients ordered in rows beginning with the splinepart for the last interval(p,p-1,..,1) and in each row with the coefficient for the highest polynomial degree (x^p, x^p-1,...,x^0)
1444  !for clamped splines in the first and last degree-1 intervals the splineparts are different
1445  case (0)
1446  self%poly_coeffs = reshape((/1._f64/), (/1, 1/))
1447  self%poly_coeffs_fp = reshape((/1._f64/), (/1, 1/))
1448  self%poly_coeffs_fpa = reshape((/1._f64, 0._f64/), (/2, 1/))
1449  case (1)
1450  self%poly_coeffs = reshape((/-1._f64, 1._f64, 1._f64, 0._f64/), (/2, 2/))
1451  self%poly_coeffs_fp = reshape((/-inv_2, 1._f64, inv_2, 0._f64/), (/2, 2/))
1452  self%poly_coeffs_fpa = reshape((/-inv_2, 1._f64, inv_2, inv_2, 0._f64, 0._f64/), (/3, 2/))
1453  self%poly_coeffs_fd = reshape((/-1._f64, 1._f64/), (/1, 2/))
1454 
1455  case (2)
1456  self%poly_coeffs = reshape((/inv_2, -1._f64, inv_2, &
1457  -1._f64, 1._f64, inv_2, &
1458  inv_2, 0._f64, 0._f64/), (/3, 3/))
1459 
1460  self%poly_coeffs_fp = reshape((/inv_6, -inv_2, inv_2, &
1461  -inv_3, inv_2, inv_2, &
1462  inv_6, 0._f64, 0._f64/), (/3, 3/))
1463 
1464  self%poly_coeffs_fpa = reshape((/inv_6, -inv_2, inv_2, 5.0_f64*inv_6, &
1465  -inv_3, inv_2, inv_2, inv_6, &
1466  inv_6, 0._f64, 0._f64, 0._f64/), (/4, 3/))
1467 
1468  self%poly_coeffs_fd = reshape((/1._f64, -1._f64, &
1469  -2._f64, 1._f64, &
1470  1._f64, 0._f64/), (/2, 3/))
1471 
1472  ! first interval
1473  self%poly_coeffs_boundary_left(:, 1, 1) = [1._f64, -2._f64, 1._f64]
1474  self%poly_coeffs_boundary_left(:, 2, 1) = [-1.5_f64, 2._f64, 0._f64]
1475  self%poly_coeffs_boundary_left(:, 3, 1) = [inv_2, 0._f64, 0._f64]
1476 
1477  self%poly_coeffs_fd_boundary_left(:, 1, 1) = [2._f64, -2._f64]
1478  self%poly_coeffs_fd_boundary_left(:, 2, 1) = [-3._f64, 2._f64]
1479  self%poly_coeffs_fd_boundary_left(:, 3, 1) = [1._f64, 0._f64]
1480 
1481  ! last interval
1482  self%poly_coeffs_boundary_right(:, 1, 1) = [0.5_f64, -1._f64, 0.5_f64]
1483  self%poly_coeffs_boundary_right(:, 2, 1) = [-1.5_f64, 1._f64, 0.5_f64]
1484  self%poly_coeffs_boundary_right(:, 3, 1) = [1.0_f64, 0._f64, 0._f64]
1485 
1486  self%poly_coeffs_fd_boundary_right(:, 1, 1) = [1._f64, -1._f64]
1487  self%poly_coeffs_fd_boundary_right(:, 2, 1) = [-3._f64, 1._f64]
1488  self%poly_coeffs_fd_boundary_right(:, 3, 1) = [2.0_f64, 0._f64]
1489 
1490  if (self%boundary_conditions == sll_p_boundary_clamped_square) then
1491  self%poly_coeffs_boundary_left(:, 1, 1) = 0._f64 !perfect conductor boundary condition
1492  self%poly_coeffs_boundary_right(:, 3, 1) = 0._f64!perfect conductor boundary condition
1493  end if
1494 
1495  case (3)
1496  self%poly_coeffs = reshape((/-inv_6, inv_2, -inv_2, inv_6, &
1497  inv_2, -1._f64, 0._f64, 4._f64*inv_6, &
1498  -inv_2, inv_2, inv_2, inv_6, &
1499  inv_6, 0._f64, 0._f64, 0._f64/), (/4, 4/))
1500 
1501  self%poly_coeffs_fp = reshape((/-inv_24, inv_6, -inv_4, inv_6, &
1502  inv_8, -inv_3, 0._f64, 4._f64*inv_6, &
1503  -inv_8, inv_6, inv_4, inv_6, &
1504  inv_24, 0._f64, 0._f64, 0._f64/), (/4, 4/))
1505 
1506  self%poly_coeffs_fpa = reshape((/-inv_24, inv_6, -inv_4, inv_6, 23._f64*inv_24, &
1507  inv_8, -inv_3, 0._f64, 4._f64*inv_6, inv_2, &
1508  -inv_8, inv_6, inv_4, inv_6, inv_24, &
1509  inv_24, 0._f64, 0._f64, 0._f64, 0._f64/), (/5, 4/))
1510 
1511  self%poly_coeffs_fd = reshape((/-inv_2, 1._f64, -inv_2, &
1512  3._f64*inv_2, -2._f64, 0._f64, &
1513  -3._f64*inv_2, 1._f64, inv_2, &
1514  inv_2, 0._f64, 0._f64/), (/3, 4/))
1515 
1516  ! first interval
1517  self%poly_coeffs_boundary_left(:, 1, 1) = [-1.0_f64, 3.0_f64, -3.0_f64, 1.0_f64]
1518  self%poly_coeffs_boundary_left(:, 2, 1) = [1.75_f64, -4.5_f64, 3.0_f64, 0._f64]
1519  self%poly_coeffs_boundary_left(:, 3, 1) = [-11._f64/12._f64, 1.5_f64, 0.0_f64, 0.0_f64]
1520  self%poly_coeffs_boundary_left(:, 4, 1) = [inv_6, 0._f64, 0._f64, 0._f64]
1521 
1522  self%poly_coeffs_fd_boundary_left(:, 1, 1) = [-3.0_f64, 6.0_f64, -3.0_f64]
1523  self%poly_coeffs_fd_boundary_left(:, 2, 1) = [5.25_f64, -9._f64, 3.0_f64]
1524  self%poly_coeffs_fd_boundary_left(:, 3, 1) = [-11._f64/4._f64, 3._f64, 0.0_f64]
1525  self%poly_coeffs_fd_boundary_left(:, 4, 1) = [inv_2, 0._f64, 0._f64]
1526 
1527  ! second interval
1528  self%poly_coeffs_boundary_left(:, 1, 2) = [-0.25_f64, 0.75_f64, -0.75_f64, 0.25_f64]
1529  self%poly_coeffs_boundary_left(:, 2, 2) = [7._f64/12._f64, -1.25_f64, 0.25_f64, 7._f64/12._f64]
1530  self%poly_coeffs_boundary_left(:, 3, 2) = [-0.5_f64, 0.5_f64, 0.5_f64, inv_6]
1531  self%poly_coeffs_boundary_left(:, 4, 2) = [inv_6, 0._f64, 0._f64, 0._f64]
1532 
1533  self%poly_coeffs_fd_boundary_left(:, 1, 2) = [-0.75_f64, 1.5_f64, -0.75_f64]
1534  self%poly_coeffs_fd_boundary_left(:, 2, 2) = [7._f64/4._f64, -2.5_f64, 0.25_f64]
1535  self%poly_coeffs_fd_boundary_left(:, 3, 2) = [-1.5_f64, 1._f64, 0.5_f64]
1536  self%poly_coeffs_fd_boundary_left(:, 4, 2) = [inv_2, 0._f64, 0._f64]
1537 
1538  ! second last interval
1539  self%poly_coeffs_boundary_right(:, 1, 1) = [-inv_6, inv_2, -inv_2, inv_6]
1540  self%poly_coeffs_boundary_right(:, 2, 1) = [inv_2, -1._f64, 0._f64, 4._f64*inv_6]
1541  self%poly_coeffs_boundary_right(:, 3, 1) = [-7.0_f64/12.0_f64, 0.5_f64, 0.5_f64, inv_6]
1542  self%poly_coeffs_boundary_right(:, 4, 1) = [0.25_f64, 0.0_f64, 0.0_f64, 0.0_f64]
1543 
1544  self%poly_coeffs_fd_boundary_right(:, 1, 1) = [-inv_2, 1._f64, -inv_2]
1545  self%poly_coeffs_fd_boundary_right(:, 2, 1) = [3._f64*inv_2, -2._f64, 0._f64]
1546  self%poly_coeffs_fd_boundary_right(:, 3, 1) = [-7.0_f64/4.0_f64, 1._f64, 0.5_f64]
1547  self%poly_coeffs_fd_boundary_right(:, 4, 1) = [0.75_f64, 0.0_f64, 0.0_f64]
1548 
1549  ! last interval
1550  self%poly_coeffs_boundary_right(:, 1, 2) = [-inv_6, inv_2, -inv_2, inv_6]
1551  self%poly_coeffs_boundary_right(:, 2, 2) = [11.0_f64/12.0_f64, -1.25_f64, -0.25_f64, 7.0_f64/12.0_f64]
1552  self%poly_coeffs_boundary_right(:, 3, 2) = [-1.75_f64, 0.75_f64, 0.75_f64, 0.25_f64]
1553  self%poly_coeffs_boundary_right(:, 4, 2) = [1.0_f64, 0.0_f64, 0.0_f64, 0.0_f64]
1554 
1555  self%poly_coeffs_fd_boundary_right(:, 1, 2) = [-inv_2, 1._f64, -inv_2]
1556  self%poly_coeffs_fd_boundary_right(:, 2, 2) = [11.0_f64/4.0_f64, -2.5_f64, -0.25_f64]
1557  self%poly_coeffs_fd_boundary_right(:, 3, 2) = [-5.25_f64, 1.5_f64, 0.75_f64]
1558  self%poly_coeffs_fd_boundary_right(:, 4, 2) = [3.0_f64, 0.0_f64, 0.0_f64]
1559 
1560  if (self%boundary_conditions == sll_p_boundary_clamped_cubic) then
1561  self%poly_coeffs_boundary_left(:, 1, 1) = 0._f64 !perfect conductor boundary condition
1562  self%poly_coeffs_boundary_right(:, 4, 2) = 0._f64 !perfect conductor boundary condition
1563  end if
1564  case (4)
1565  self%poly_coeffs = reshape((/inv_24, -inv_6, inv_4, -inv_6, inv_24 &
1566  , -inv_6, inv_2, -inv_4, -inv_2, 11._f64*inv_24 &
1567  , inv_4, -inv_2, -inv_4, inv_2, 11._f64*inv_24 &
1568  , -inv_6, inv_6, inv_4, inv_6, inv_24 &
1569  , inv_24, 0._f64, 0._f64, 0._f64, 0._f64/), (/5, 5/))
1570 
1571  self%poly_coeffs_fp = reshape((/inv_120, -inv_24, inv_12, -inv_12, inv_24 &
1572  , -inv_30, inv_8, -inv_12, -inv_4, 11._f64*inv_24 &
1573  , inv_20, -inv_8, -inv_12, inv_4, 11._f64*inv_24 &
1574  , -inv_30, inv_24, inv_12, inv_12, inv_24 &
1575  , inv_120, 0._f64, 0._f64, 0._f64, 0._f64/), (/5, 5/))
1576 
1577  self%poly_coeffs_fpa = reshape((/inv_120, -inv_24, inv_12, -inv_12, inv_24, 119._f64*inv_120 &
1578  , -inv_30, inv_8, -inv_12, -inv_4, 11._f64*inv_24, 93._f64*inv_120 &
1579  , inv_20, -inv_8, -inv_12, inv_4, 11._f64*inv_24, 27._f64*inv_120 &
1580  , -inv_30, inv_24, inv_12, inv_12, inv_24, inv_120 &
1581  , inv_120, 0._f64, 0._f64, 0._f64, 0._f64, 0.0_f64/), (/6, 5/))
1582 
1583  case (5)
1584  self%poly_coeffs = reshape((/-inv_120, inv_24, -inv_12, inv_12, -inv_24, inv_120 &
1585  , inv_24, -inv_6, inv_6, inv_6, -5._f64*inv_12, 26._f64*inv_120 &
1586  , -inv_12, inv_4, 0._f64, -inv_2, 0._f64, 11._f64*inv_20 &
1587  , inv_12, -inv_6, -inv_6, inv_6, 5._f64*inv_12, 26._f64*inv_120 &
1588  , -inv_24, inv_24, inv_12, inv_12, inv_24, inv_120 &
1589  , inv_120, 0._f64, 0._f64, 0._f64, 0._f64, 0._f64/), (/6, 6/))
1590 
1591  self%poly_coeffs_fp = reshape((/-inv_720, inv_120, -inv_48, inv_36, -inv_48, inv_120 &
1592  , inv_144, -inv_30, inv_24, inv_18, -5._f64*inv_24, 26._f64*inv_120 &
1593  , -inv_72, inv_20, 0._f64, -inv_6, 0._f64, 11._f64*inv_20 &
1594  , inv_72, -inv_30, -inv_24, inv_18, 5._f64*inv_24, 26._f64*inv_120 &
1595  , -inv_144, inv_120, inv_48, inv_36, inv_48, inv_120 &
1596  , inv_720, 0._f64, 0._f64, 0._f64, 0._f64, 0._f64/), (/6, 6/))
1597 
1598  self%poly_coeffs_fpa = reshape((/-inv_720, inv_120, -inv_48, inv_36, -inv_48, inv_120, 719._f64*inv_720 &
1599  , inv_144, -inv_30, inv_24, inv_18, -5._f64*inv_24, 26._f64*inv_120, 662._f64*inv_720 &
1600  , -inv_72, inv_20, 0._f64, -inv_6, 0._f64, 11._f64*inv_20, inv_2 &
1601  , inv_72, -inv_30, -inv_24, inv_18, 5._f64*inv_24, 26._f64*inv_120, 58._f64*inv_720 &
1602  , -inv_144, inv_120, inv_48, inv_36, inv_48, inv_120, inv_720 &
1603  , inv_720, 0._f64, 0._f64, 0._f64, 0._f64, 0._f64, 0._f64/), (/7, 6/))
1604  case (6)
1605  self%poly_coeffs = reshape((/inv_720, -inv_120, inv_48, -inv_36, inv_48, -inv_120, inv_720 &
1606  , -inv_120, inv_24, -inv_16, -inv_36, 135._f64*inv_720, -150._f64*inv_720, 57._f64*inv_720 &
1607  , inv_48, -inv_12, inv_24, 160._f64*inv_720, -150._f64*inv_720, -inv_3, 302._f64*inv_720 &
1608  , -inv_36, inv_12, inv_24, -160._f64*inv_720, -150._f64*inv_720, inv_3, 302._f64*inv_720 &
1609  , inv_48, -inv_24, -inv_16, inv_36, 135._f64*inv_720, 150._f64*inv_720, 57._f64*inv_720 &
1610  , -inv_120, inv_120, inv_48, inv_36, inv_48, inv_120, inv_720 &
1611  , inv_720, 0._f64, 0._f64, 0._f64, 0._f64, 0._f64, 0._f64/), (/7, 7/))
1612 
1613  self%poly_coeffs_fp = reshape((/inv_5040, -inv_720, inv_240, -inv_144, inv_144, -inv_240, inv_720 &
1614  , -6._f64*inv_5040, inv_144, -inv_80, -inv_144, inv_16, -75._f64*inv_720, 57._f64*inv_720 &
1615  , 15._f64*inv_5040, -inv_72, inv_120, inv_18, -50._f64*inv_720, -inv_6, 302._f64*inv_720 &
1616  , -20._f64*inv_5040, inv_72, inv_120, -inv_18, -50._f64*inv_720, inv_6, 302._f64*inv_720 &
1617  , 15._f64*inv_5040, -inv_144, -inv_80, inv_144, inv_16, 75._f64*inv_720, 57._f64*inv_720 &
1618  , -6.0_f64*inv_5040, inv_720, inv_240, inv_144, inv_144, inv_240, inv_720 &
1619  , inv_5040, 0._f64, 0._f64, 0._f64, 0._f64, 0._f64, 0.0_f64/), (/7, 7/))
1620 
1621  self%poly_coeffs_fpa = reshape((/inv_5040, -inv_720, inv_240, -inv_144, inv_144, -inv_240, inv_720, 5039._f64*inv_5040 &
1622  , -6._f64*inv_5040, inv_144, -inv_80, -inv_144, inv_16, -75._f64*inv_720, 57._f64*inv_720, 4919._f64*inv_5040 &
1623  , 15._f64*inv_5040, -inv_72, inv_120, inv_18, -50._f64*inv_720, -inv_6, 302._f64*inv_720, 3728._f64*inv_5040 &
1624  , -20._f64*inv_5040, inv_72, inv_120, -inv_18, -50._f64*inv_720, inv_6, 302._f64*inv_720, 1312._f64*inv_5040 &
1625  , 15._f64*inv_5040, -inv_144, -inv_80, inv_144, inv_16, 75._f64*inv_720, 57._f64*inv_720, 121._f64*inv_5040 &
1626  , -6.0_f64*inv_5040, inv_720, inv_240, inv_144, inv_144, inv_240, inv_720, 1._f64*inv_5040 &
1627  , inv_5040, 0._f64, 0._f64, 0._f64, 0._f64, 0._f64, 0.0_f64, 0._f64/), (/8, 7/))
1628  case default
1629  sll_error('sll_s_spline_pp_init', 'Not implemented')
1630  end select
1631 
1632  if (self%boundary_conditions > 0 .and. degree > 3) then
1633  sll_error('sll_s_spline_pp_init_1d', 'Specified boundary conditions not implemented for given degree.')
1634  end if
1635 
1636  end subroutine sll_s_spline_pp_init_1d
1637 
1639  subroutine sll_s_spline_pp_init_2d(self, degree, n_cells, boundary)
1640  type(sll_t_spline_pp_2d) self
1641  sll_int32, intent(in) :: degree(2)
1642  sll_int32, intent(in) :: n_cells(2)
1643  sll_int32, intent(in), optional :: boundary(2)
1644 
1645  if (present(boundary)) then
1646  call sll_s_spline_pp_init_1d(self%spline1, degree(1), n_cells(1), boundary(1))
1647  call sll_s_spline_pp_init_1d(self%spline2, degree(2), n_cells(2), boundary(2))
1648  else
1649  call sll_s_spline_pp_init_1d(self%spline1, degree(1), n_cells(1))
1650  call sll_s_spline_pp_init_1d(self%spline2, degree(2), n_cells(2))
1651  end if
1652 
1653  end subroutine sll_s_spline_pp_init_2d
1654 
1656  subroutine sll_s_spline_pp_init_3d(self, degree, n_cells, boundary)
1657  type(sll_t_spline_pp_3d) self
1658  sll_int32, intent(in) :: degree(3)
1659  sll_int32, intent(in) :: n_cells(3)
1660  sll_int32, intent(in), optional :: boundary(3)
1661  sll_int32 :: ierr
1662 
1663  if (present(boundary)) then
1664  call sll_s_spline_pp_init_1d(self%spline1, degree(1), n_cells(1), boundary(1))
1665  call sll_s_spline_pp_init_1d(self%spline2, degree(2), n_cells(2), boundary(2))
1666  call sll_s_spline_pp_init_1d(self%spline3, degree(3), n_cells(3), boundary(3))
1667  else
1668  call sll_s_spline_pp_init_1d(self%spline1, degree(1), n_cells(1))
1669  call sll_s_spline_pp_init_1d(self%spline2, degree(2), n_cells(2))
1670  call sll_s_spline_pp_init_1d(self%spline3, degree(3), n_cells(3))
1671  end if
1672  sll_allocate(self%scratch_b((degree(2) + 1)*(degree(3) + 1)), ierr)
1673  sll_assert(ierr == 0)
1674  sll_allocate(self%scratch_p((degree(2) + 1)*(degree(3) + 1)), ierr)
1675  sll_assert(ierr == 0)
1676  sll_allocate(self%scratch_pp((degree(1) + 1)*(degree(2) + 1)*(degree(3) + 1)), ierr)
1677  sll_assert(ierr == 0)
1678  sll_allocate(self%scratch_coeffs((degree(1) + 1)*(degree(2) + 1)*(degree(3) + 1), n_cells(1)*n_cells(2)*n_cells(3)), ierr)
1679  sll_assert(ierr == 0)
1680 
1681  end subroutine sll_s_spline_pp_init_3d
1682 
1684  subroutine sll_s_spline_pp_free_1d(self)
1685  type(sll_t_spline_pp_1d), intent(inout) :: self
1686 
1687  deallocate (self%poly_coeffs)
1688  deallocate (self%poly_coeffs_fp)
1689  deallocate (self%scratch_b)
1690  deallocate (self%scratch_p)
1691 
1692  end subroutine sll_s_spline_pp_free_1d
1693 
1695  subroutine sll_s_spline_pp_free_2d(self)
1696  type(sll_t_spline_pp_2d), intent(inout) :: self
1697 
1698  call sll_s_spline_pp_free_1d(self%spline1)
1699  call sll_s_spline_pp_free_1d(self%spline2)
1700 
1701  end subroutine sll_s_spline_pp_free_2d
1702 
1704  subroutine sll_s_spline_pp_free_3d(self)
1705  type(sll_t_spline_pp_3d), intent(inout) :: self
1706 
1707  call sll_s_spline_pp_free_1d(self%spline1)
1708  call sll_s_spline_pp_free_1d(self%spline2)
1709  call sll_s_spline_pp_free_1d(self%spline3)
1710  deallocate (self%scratch_b)
1711  deallocate (self%scratch_p)
1712  deallocate (self%scratch_pp)
1713 
1714  end subroutine sll_s_spline_pp_free_3d
1715 
1716 end module sll_m_splines_pp
sll_m_splines_pp
Splines in pp form.
Definition: sll_m_splines_pp.F90:8
sll_m_splines_pp::inv_72
real(kind=f64), parameter inv_72
Definition: sll_m_splines_pp.F90:84
sll_m_splines_pp::inv_12
real(kind=f64), parameter inv_12
Definition: sll_m_splines_pp.F90:75
sll_m_splines_pp::sll_s_spline_pp_b_to_pp_3d_clamped_2full
subroutine, public sll_s_spline_pp_b_to_pp_3d_clamped_2full(self, n_cells, b_coeffs, pp_coeffs)
Convert 3d spline in B form to spline in pp form for clamped spline.
Definition: sll_m_splines_pp.F90:975
sll_m_splines_pp::sll_s_spline_pp_init_3d
subroutine, public sll_s_spline_pp_init_3d(self, degree, n_cells, boundary)
Initialize sll_t_spline_pp_3d object.
Definition: sll_m_splines_pp.F90:1657
sll_m_splines_pp::inv_6
real(kind=f64), parameter inv_6
Definition: sll_m_splines_pp.F90:72
sll_m_splines_pp::sll_f_spline_pp_horner_derivative_3d
real(kind=f64) function, public sll_f_spline_pp_horner_derivative_3d(degree, pp_coeffs, x, indices, n_cells, component)
Perform a 3d hornerschema on the pp_coeffs at the indices.
Definition: sll_m_splines_pp.F90:1347
sll_m_splines_pp::sll_f_spline_pp_horner_2d
real(kind=f64) function, public sll_f_spline_pp_horner_2d(degree, pp_coeffs, x, indices, n_cells)
Perform a 2d hornerschema on the pp_coeffs at the indices.
Definition: sll_m_splines_pp.F90:1209
sll_m_splines_pp::sll_s_spline_pp_b_to_pp_3d_clamped
subroutine, public sll_s_spline_pp_b_to_pp_3d_clamped(self, n_cells, b_coeffs, pp_coeffs)
Convert 3d spline in B form to spline in pp form for clamped spline.
Definition: sll_m_splines_pp.F90:958
sll_m_splines_pp::inv_20
real(kind=f64), parameter inv_20
Definition: sll_m_splines_pp.F90:78
sll_m_splines_pp::sll_s_spline_pp_b_to_pp_1d_cell
subroutine sll_s_spline_pp_b_to_pp_1d_cell(self, b_coeffs, pp_coeffs)
Convert 1d spline in B form in a cell to spline in pp form with periodic boundary conditions.
Definition: sll_m_splines_pp.F90:228
sll_m_splines_pp::sll_s_spline_pp_b_to_pp_3d_cella2f
subroutine sll_s_spline_pp_b_to_pp_3d_cella2f(self, n_cells, b_coeffs, pp_coeffs, i, j, k)
Convert 3d spline in B form in a cell to spline in pp form with clamped boundary in all three directi...
Definition: sll_m_splines_pp.F90:725
sll_m_splines_pp::sll_s_spline_pp_horner_m_3d
subroutine, public sll_s_spline_pp_horner_m_3d(self, val, degree, x)
Perform three times a 1d hornerschema on the poly_coeffs.
Definition: sll_m_splines_pp.F90:1143
sll_m_splines_pp::inv_30
real(kind=f64), parameter inv_30
Definition: sll_m_splines_pp.F90:80
sll_m_splines_pp::inv_720
real(kind=f64), parameter inv_720
Definition: sll_m_splines_pp.F90:90
sll_m_splines_pp::sll_s_spline_evaluate_basis_b_form_3d_clamped
subroutine, public sll_s_spline_evaluate_basis_b_form_3d_clamped(self, n_cells, b_coeffs, val)
Definition: sll_m_splines_pp.F90:1054
sll_m_splines_pp::sll_p_boundary_clamped
integer(kind=i32), parameter, public sll_p_boundary_clamped
Parameter specifying clamped boundary conditions.
Definition: sll_m_splines_pp.F90:62
sll_m_splines_pp::inv_24
real(kind=f64), parameter inv_24
Definition: sll_m_splines_pp.F90:79
sll_m_splines_pp::sll_s_spline_pp_b_to_pp_1d_clamped_clampeddiri
subroutine sll_s_spline_pp_b_to_pp_1d_clamped_clampeddiri(self, n_cells, b_coeffs, pp_coeffs)
Convert 1d spline in B form to spline in pp form for clamped spline and Dirichlet conditions on the r...
Definition: sll_m_splines_pp.F90:182
sll_m_splines_pp::sll_s_spline_pp_b_to_pp_3d
subroutine, public sll_s_spline_pp_b_to_pp_3d(self, n_cells, b_coeffs, pp_coeffs)
Convert 3d spline in B form to spline in pp form with periodic boundary conditions.
Definition: sll_m_splines_pp.F90:529
sll_m_splines_pp::sll_s_spline_pp_b_to_pp_1d_cella
subroutine sll_s_spline_pp_b_to_pp_1d_cella(degree, pp_b, b_coeffs, pp_coeffs)
Convert 1d spline in B form in a cell to spline in pp form for specified pp_coefficients of the b_spl...
Definition: sll_m_splines_pp.F90:246
sll_m_splines_pp::sll_f_spline_pp_horner_3d
real(kind=f64) function, public sll_f_spline_pp_horner_3d(degree, pp_coeffs, x, indices, n_cells)
Perform a 3d hornerschema on the pp_coeffs at the indices.
Definition: sll_m_splines_pp.F90:1227
sll_m_splines_pp::sll_s_spline_pp_horner_primitive_1d
subroutine, public sll_s_spline_pp_horner_primitive_1d(val, degree, pp_coeffs, x)
Perform a 1d hornerschema on the pp_coeffs evaluate at x.
Definition: sll_m_splines_pp.F90:1162
sll_m_splines_pp::inv_120
real(kind=f64), parameter inv_120
Definition: sll_m_splines_pp.F90:86
sll_m_splines_pp::sll_s_spline_pp_horner_m_2d
subroutine, public sll_s_spline_pp_horner_m_2d(self, val, degree, x)
Perform two times a 1d hornerschema on the poly_coeffs.
Definition: sll_m_splines_pp.F90:1127
sll_m_splines_pp::sll_s_spline_pp_free_3d
subroutine, public sll_s_spline_pp_free_3d(self)
Destructor 3d.
Definition: sll_m_splines_pp.F90:1705
sll_m_splines_pp::inv_48
real(kind=f64), parameter inv_48
Definition: sll_m_splines_pp.F90:82
sll_m_splines_pp::inv_16
real(kind=f64), parameter inv_16
Definition: sll_m_splines_pp.F90:76
sll_m_splines_pp::sll_s_spline_evaluate_basis_b_form_3d_periodic
subroutine, public sll_s_spline_evaluate_basis_b_form_3d_periodic(self, n_cells, b_coeffs, val)
Definition: sll_m_splines_pp.F90:1066
sll_m_splines_pp::sll_s_spline_pp_b_to_pp_3d_cella
subroutine sll_s_spline_pp_b_to_pp_3d_cella(self, n_cells, b_coeffs, pp_coeffs, i, j, k)
Convert 3d spline in B form in a cell to spline in pp form with clamped boundary in first direction a...
Definition: sll_m_splines_pp.F90:631
sll_m_splines_pp::inv_80
real(kind=f64), parameter inv_80
Definition: sll_m_splines_pp.F90:85
sll_m_splines_pp::inv_144
real(kind=f64), parameter inv_144
Definition: sll_m_splines_pp.F90:88
sll_m_splines_pp::sll_s_spline_pp_b_to_pp_2d_periodic
subroutine, public sll_s_spline_pp_b_to_pp_2d_periodic(self, n_cells, b_coeffs, pp_coeffs)
Convert 2d spline in B form to spline in pp form with periodic boundary conditions This is a special ...
Definition: sll_m_splines_pp.F90:266
sll_m_splines_pp::sll_s_spline_pp_init_1d
subroutine, public sll_s_spline_pp_init_1d(self, degree, n_cells, boundary)
Initialize sll_t_spline_pp_1d object (Set poly_coeffs depending on degree)
Definition: sll_m_splines_pp.F90:1387
sll_m_splines_pp::sll_s_spline_pp_free_1d
subroutine, public sll_s_spline_pp_free_1d(self)
Destructor 1d.
Definition: sll_m_splines_pp.F90:1685
sll_m_splines_pp::sll_s_spline_pp_b_to_pp_3d_cell
subroutine sll_s_spline_pp_b_to_pp_3d_cell(self, n_cells, b_coeffs, pp_coeffs, i, j, k)
Convert 3d spline in B form in a cell to spline in pp form with periodic boundary conditions.
Definition: sll_m_splines_pp.F90:547
sll_m_splines_pp::inv_5040
real(kind=f64), parameter inv_5040
Definition: sll_m_splines_pp.F90:91
sll_m_splines_pp::sll_s_spline_pp_b_to_pp_1d
subroutine, public sll_s_spline_pp_b_to_pp_1d(self, n_cells, b_coeffs, pp_coeffs)
Convert 1d spline in B form to spline in pp form with periodic boundary conditions.
Definition: sll_m_splines_pp.F90:210
sll_m_splines_pp::inv_36
real(kind=f64), parameter inv_36
Definition: sll_m_splines_pp.F90:81
sll_m_splines_pp::sll_s_spline_pp_horner_m_1d
subroutine, public sll_s_spline_pp_horner_m_1d(self, val, degree, x)
Perform a 1d hornerschema on the poly_coeffs.
Definition: sll_m_splines_pp.F90:1114
sll_m_splines_pp::sll_s_spline_pp_b_to_pp_3d_cellaf
subroutine sll_s_spline_pp_b_to_pp_3d_cellaf(self, n_cells, b_coeffs, pp_coeffs, i, j, k)
Convert 3d spline in B form in a cell to spline in pp form with clamped boundary in all three directi...
Definition: sll_m_splines_pp.F90:834
sll_m_splines_pp::inv_7
real(kind=f64), parameter inv_7
Definition: sll_m_splines_pp.F90:73
sll_m_splines_pp::sll_p_boundary_periodic
integer(kind=i32), parameter, public sll_p_boundary_periodic
Parameter to specify boundary conditions.
Definition: sll_m_splines_pp.F90:61
sll_m_splines_pp::sll_p_boundary_clamped_clampeddiri
integer(kind=i32), parameter, public sll_p_boundary_clamped_clampeddiri
Parameter specifying clamped boundary conditions with the last point left out (for Dirichlet conditio...
Definition: sll_m_splines_pp.F90:63
sll_m_splines_pp::sll_s_spline_evaluate_basis_pp_form_1d
subroutine sll_s_spline_evaluate_basis_pp_form_1d(self, n_cells, pp_coeffs, val)
Definition: sll_m_splines_pp.F90:1038
sll_m_splines_pp::sll_f_spline_pp_horner_3d_d1
real(kind=f64) function, public sll_f_spline_pp_horner_3d_d1(degree, pp_coeffs, x, indices, n_cells)
Perform a 3d hornerschema on the pp_coeffs at the indices.
Definition: sll_m_splines_pp.F90:1251
sll_m_splines_pp::inv_2
real(kind=f64), parameter inv_2
Definition: sll_m_splines_pp.F90:69
sll_m_splines_pp::sll_s_spline_pp_b_to_pp_2d_cell
subroutine sll_s_spline_pp_b_to_pp_2d_cell(spline1, spline2, n_cells, b_coeffs, pp_coeffs, i, j)
Convert 2d spline in B form in a cell to spline in pp form with periodic boundary conditions.
Definition: sll_m_splines_pp.F90:478
sll_m_splines_pp::sll_f_spline_pp_horner_1d
real(kind=f64) function, public sll_f_spline_pp_horner_1d(degree, pp_coeffs, x, index)
Perform a 1d hornerschema on the pp_coeffs at index.
Definition: sll_m_splines_pp.F90:1175
sll_m_splines_pp::inv_60
real(kind=f64), parameter inv_60
Definition: sll_m_splines_pp.F90:83
sll_m_splines_pp::sll_s_spline_pp_b_to_pp_1d_clamped
subroutine sll_s_spline_pp_b_to_pp_1d_clamped(self, n_cells, b_coeffs, pp_coeffs)
Convert 1d spline in B form to spline in pp form for clamped spline.
Definition: sll_m_splines_pp.F90:156
sll_m_splines_pp::sll_s_spline_pp_init_2d
subroutine, public sll_s_spline_pp_init_2d(self, degree, n_cells, boundary)
Initialize sll_t_spline_pp_2d object.
Definition: sll_m_splines_pp.F90:1640
sll_m_splines_pp::sll_p_boundary_clampeddiri
integer(kind=i32), parameter, public sll_p_boundary_clampeddiri
Parameter specifying clamped boundary conditions with the first and last point left out (for Dirichle...
Definition: sll_m_splines_pp.F90:64
sll_m_splines_pp::sll_s_spline_pp_b_to_pp_3d_clamped_full
subroutine sll_s_spline_pp_b_to_pp_3d_clamped_full(self, n_cells, b_coeffs, pp_coeffs)
Convert 3d spline in B form to spline in pp form for clamped spline.
Definition: sll_m_splines_pp.F90:993
sll_m_splines_pp::sll_f_spline_pp_horner_derivative_1d
real(kind=f64) function sll_f_spline_pp_horner_derivative_1d(degree, pp_coeffs, x, index)
Perform a 1d hornerschema on the pp_coeffs at index.
Definition: sll_m_splines_pp.F90:1192
sll_m_splines_pp::sll_p_boundary_clamped_square
integer(kind=i32), parameter, public sll_p_boundary_clamped_square
Parameter specifying clamped boundary conditions with square spline for the 0-form,...
Definition: sll_m_splines_pp.F90:66
sll_m_splines_pp::sll_s_spline_pp_pp_to_b_1d
subroutine, public sll_s_spline_pp_pp_to_b_1d(self, n_cells, pp_coeffs, b_coeffs)
Compute b_coeffs from coefficients of the monomials (in pp_coeffs)
Definition: sll_m_splines_pp.F90:136
sll_m_splines_pp::inv_3
real(kind=f64), parameter inv_3
Definition: sll_m_splines_pp.F90:70
sll_m_splines_pp::sll_s_spline_pp_b_to_pp_2d
subroutine, public sll_s_spline_pp_b_to_pp_2d(self, n_cells, b_coeffs, pp_coeffs)
Convert 2d spline in B form to spline in pp form.
Definition: sll_m_splines_pp.F90:285
sll_m_splines_pp::sll_s_spline_evaluate_basis_b_form_1d_periodic
subroutine, public sll_s_spline_evaluate_basis_b_form_1d_periodic(self, n_cells, b_coeffs, val)
Definition: sll_m_splines_pp.F90:1024
sll_m_splines_pp::inv_18
real(kind=f64), parameter inv_18
Definition: sll_m_splines_pp.F90:77
sll_m_splines_pp::sll_s_spline_evaluate_basis_b_form_1d_clamped
subroutine, public sll_s_spline_evaluate_basis_b_form_1d_clamped(self, n_cells, b_coeffs, val)
Definition: sll_m_splines_pp.F90:1010
sll_m_splines_pp::sll_f_spline_pp_horner_3d_d3
real(kind=f64) function, public sll_f_spline_pp_horner_3d_d3(degree, pp_coeffs, x, indices, n_cells)
Perform a 3d hornerschema on the pp_coeffs at the indices.
Definition: sll_m_splines_pp.F90:1315
sll_m_splines_pp::sll_p_boundary_clampeddiri_clamped
integer(kind=i32), parameter, public sll_p_boundary_clampeddiri_clamped
Parameter specifying clamped boundary conditions with the first point left out (for Dirichlet conditi...
Definition: sll_m_splines_pp.F90:65
sll_m_splines_pp::sll_p_boundary_clamped_cubic
integer(kind=i32), parameter, public sll_p_boundary_clamped_cubic
Parameter specifying clamped boundary conditions with a cubic spline for the 0-form,...
Definition: sll_m_splines_pp.F90:67
sll_m_splines_pp::inv_8
real(kind=f64), parameter inv_8
Definition: sll_m_splines_pp.F90:74
sll_m_splines_pp::inv_240
real(kind=f64), parameter inv_240
Definition: sll_m_splines_pp.F90:87
sll_m_splines_pp::sll_s_spline_pp_free_2d
subroutine, public sll_s_spline_pp_free_2d(self)
Destructor 2d.
Definition: sll_m_splines_pp.F90:1696
sll_m_splines_pp::inv_4
real(kind=f64), parameter inv_4
Definition: sll_m_splines_pp.F90:71
sll_m_splines_pp::sll_f_spline_pp_horner_3d_d2
real(kind=f64) function, public sll_f_spline_pp_horner_3d_d2(degree, pp_coeffs, x, indices, n_cells)
Perform a 3d hornerschema on the pp_coeffs at the indices.
Definition: sll_m_splines_pp.F90:1283
sll_m_splines_pp::sll_s_spline_evaluate_basis_pp_form_3d
subroutine sll_s_spline_evaluate_basis_pp_form_3d(self, n_cells, pp_coeffs, val)
Definition: sll_m_splines_pp.F90:1078
sll_m_splines_pp::sll_t_spline_pp_1d
arbitrary degree 1d spline
Definition: sll_m_splines_pp.F90:93
sll_m_splines_pp::sll_t_spline_pp_2d
arbitrary degree 2d spline
Definition: sll_m_splines_pp.F90:113
sll_m_splines_pp::sll_t_spline_pp_3d
arbitrary degree 3d spline
Definition: sll_m_splines_pp.F90:120
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:41 for Semi-Lagrangian Library by   doxygen 1.9.1