Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_cubic_non_uniform_splines.F90
Go to the documentation of this file.
1 
5 
6 module sll_m_cubic_non_uniform_splines
7 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8 #include "sll_assert.h"
9 #include "sll_memory.h"
10 #include "sll_working_precision.h"
11 
12  use sll_m_boundary_condition_descriptors, only: &
13  sll_p_hermite, &
14  sll_p_periodic
15 
16  use sll_m_tridiagonal, only: &
17  sll_s_setup_cyclic_tridiag, &
18  sll_o_solve_cyclic_tridiag
19 
20  implicit none
21 
22  public :: &
23  sll_s_compute_spline_nonunif, &
24  sll_s_compute_spline_nonunif_1d_periodic_aux2, &
25  sll_t_cubic_nonunif_spline_1d, &
26  sll_s_interpolate_array_value_nonunif, &
27  sll_s_interpolate_array_value_nonunif_aux, &
28  sll_f_new_cubic_nonunif_spline_1d, &
29  sll_s_setup_spline_nonunif_1d_periodic_aux, &
30  sll_o_delete
31 
32  private
33 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
34 
36  type :: sll_t_cubic_nonunif_spline_1d
37  sll_int32 :: n_cells
38  sll_real64, dimension(:), pointer :: node_positions
40  sll_real64, dimension(:), pointer :: buf
41  sll_int32 :: size_buf
42  sll_int32, dimension(:), pointer :: ibuf
43  sll_int32 :: size_ibuf
44  sll_real64, dimension(:), pointer :: coeffs
45  sll_real64 :: xmin
46  sll_real64 :: xmax
47  sll_int32 :: bc_type
48  LOGICAL :: is_setup
49  sll_real64 :: slope_l
50  sll_real64 :: slope_r
51  contains
52 
53  procedure :: init => sll_s_cubic_nonunif_spline_1d_init
54 
55  end type sll_t_cubic_nonunif_spline_1d
56 
57  interface sll_o_delete
58  module procedure delete_cubic_nonunif_spline_1d
59  end interface sll_o_delete
60 
61 contains ! ****************************************************************
62 
64  function sll_f_new_cubic_nonunif_spline_1d(n_cells, bc_type)
65 
66  type(sll_t_cubic_nonunif_spline_1d), pointer :: sll_f_new_cubic_nonunif_spline_1d
67  sll_int32, intent(in) :: n_cells
68  sll_int32, intent(in) :: bc_type
69  sll_int32 :: ierr
70 
71  sll_allocate(sll_f_new_cubic_nonunif_spline_1d, ierr)
72  call sll_f_new_cubic_nonunif_spline_1d%init(n_cells, bc_type)
73 
74  end function sll_f_new_cubic_nonunif_spline_1d
75 
77  subroutine sll_s_cubic_nonunif_spline_1d_init(self, n_cells, bc_type)
78 
79  class(sll_t_cubic_nonunif_spline_1d) :: self
80  sll_int32, intent(in) :: n_cells
81  sll_int32, intent(in) :: bc_type
82  sll_int32 :: ierr, size_buf, size_ibuf
83 
84  self%n_cells = n_cells
85  self%bc_type = bc_type
86 
87  sll_allocate(self%node_positions(-2:n_cells + 2), ierr)
88  size_buf = 10*(n_cells + 1)
89  size_ibuf = n_cells + 1
90  self%size_buf = size_buf
91  sll_allocate(self%buf(size_buf), ierr)
92  self%size_ibuf = size_ibuf
93  sll_allocate(self%ibuf(size_ibuf), ierr)
94  sll_allocate(self%coeffs(-1:n_cells + 1), ierr)
95  self%is_setup = .false.
96  self%slope_L = 0.0_f64
97  self%slope_R = 0.0_f64
98 
99  end subroutine sll_s_cubic_nonunif_spline_1d_init
100 
102  subroutine delete_cubic_nonunif_spline_1d(spline, ierr)
103  type(sll_t_cubic_nonunif_spline_1d), pointer :: spline
104  sll_int32 :: ierr
105  sll_assert(associated(spline))
106  sll_deallocate(spline%node_positions, ierr)
107  sll_deallocate(spline%buf, ierr)
108  sll_deallocate(spline%coeffs, ierr)
109  sll_deallocate(spline%ibuf, ierr)
110  sll_deallocate(spline, ierr)
111  end subroutine delete_cubic_nonunif_spline_1d
112 
114  subroutine sll_s_compute_spline_nonunif(f, spline, node_positions, sl, sr)
115  sll_real64, dimension(:), intent(in) :: f ! data to be fit
116  type(sll_t_cubic_nonunif_spline_1d), pointer :: spline
117  sll_real64, dimension(:), intent(in), optional :: node_positions ! non uniform mesh
118  sll_real64, intent(in), optional :: sl
119  sll_real64, intent(in), optional :: sr
120  sll_int32 ::n_cells, bc_type
121  sll_real64, dimension(:), pointer :: node_pos
122  sll_real64 :: length
123 
124  !check that the spline is initialized
125  sll_assert(associated(spline))
126  n_cells = spline%n_cells
127  bc_type = spline%bc_type
128  !copy on non_uniform mesh in spline object
129 
130  if (present(node_positions)) then
131  node_pos => spline%node_positions
132  spline%xmin = node_positions(1)
133  spline%xmax = node_positions(n_cells + 1)
134  length = spline%xmax - spline%xmin
135  node_pos(1:n_cells - 1) = (node_positions(2:n_cells) - node_positions(1))/length
136  node_pos(0) = 0.0_f64
137  node_pos(n_cells) = 1.0_f64
138  select case (bc_type)
139  case (sll_p_periodic)
140  call sll_s_setup_spline_nonunif_1d_periodic_aux(node_pos, n_cells, spline%buf, spline%ibuf)
141  case (sll_p_hermite)
142  call setup_spline_nonunif_1d_hermite_aux(node_pos, n_cells, spline%buf, spline%ibuf)
143  case default
144  print *, 'ERROR: compute_spline_nonunif_1D(): not recognized boundary condition'
145  stop
146  end select
147  spline%is_setup = .true.
148  end if
149 
150  if (.not. (spline%is_setup)) then
151  ! FIXME: THROW ERROR
152  print *, 'ERROR: compute_spline_nonunif_1D_periodic(): '
153  print *, 'node_positions are to be given first'
154  stop
155  end if
156 
157  select case (bc_type)
158  case (sll_p_periodic)
159  call compute_spline_nonunif_1d_periodic(f, spline)
160  case (sll_p_hermite)
161  if (present(sl)) then
162  spline%slope_L = sl
163  end if
164  if (present(sr)) then
165  spline%slope_R = sr
166  end if
167  call compute_spline_nonunif_1d_hermite(f, spline)
168  case default
169  print *, 'ERROR: compute_spline_nonunif_1D(): not recognized boundary condition'
170  stop
171  end select
172  end subroutine sll_s_compute_spline_nonunif
173 
175  subroutine compute_spline_nonunif_1d_periodic(f, spline)
176  sll_real64, dimension(:), intent(in), target :: f ! data to be fit
177  type(sll_t_cubic_nonunif_spline_1d), pointer :: spline
178  sll_real64, dimension(:), pointer :: buf, fp, coeffs
179  sll_int32, dimension(:), pointer :: ibuf
180  sll_int32 :: nc
181  if (.not. associated(spline)) then
182  ! FIXME: THROW ERROR
183  print *, 'ERROR: compute_spline_nonunif_1D_periodic(): ', &
184  'uninitialized spline object passed as argument. Exiting... '
185  stop
186  end if
187  if (.not. (size(f) .ge. spline%n_cells + 1)) then
188  ! FIXME: THROW ERROR
189  print *, 'ERROR: compute_spline_nonunif_1D_periodic(): '
190  write (*, '(a, i8, a, i8)') 'spline object needs data of size >= ', &
191  spline%n_cells + 1, ' . Passed size: ', size(f)
192  stop
193  end if
194  nc = spline%n_cells
195 
196  fp => f
197  buf => spline%buf
198  ibuf => spline%ibuf
199  coeffs => spline%coeffs
200  call compute_spline_nonunif_1d_periodic_aux(fp, nc, buf, ibuf, coeffs)
201  end subroutine compute_spline_nonunif_1d_periodic
202 
204  subroutine compute_spline_nonunif_1d_hermite(f, spline)
205  sll_real64, dimension(:), intent(in), target :: f ! data to be fit
206  type(sll_t_cubic_nonunif_spline_1d), pointer :: spline
207  sll_real64, dimension(:), pointer :: buf, fp, coeffs
208  sll_int32, dimension(:), pointer :: ibuf
209  sll_int32 :: nc
210  sll_real64 :: lift(4, 2)
211  if (.not. associated(spline)) then
212  ! FIXME: THROW ERROR
213  print *, 'ERROR: compute_spline_nonunif_1D_hermite(): ', &
214  'uninitialized spline object passed as argument. Exiting... '
215  stop
216  end if
217  if (.not. (size(f) .ge. spline%n_cells + 1)) then
218  ! FIXME: THROW ERROR
219  print *, 'ERROR: compute_spline_nonunif_1D_hermite(): '
220  write (*, '(a, i8, a, i8)') 'spline object needs data of size >= ', &
221  spline%n_cells + 1, ' . Passed size: ', size(f)
222  stop
223  end if
224  nc = spline%n_cells
225 
226  fp => f
227  buf => spline%buf
228  ibuf => spline%ibuf
229  coeffs => spline%coeffs
230 
231  !lift(1,1) = 1._f64/3._f64*(spline%node_positions(1)-spline%node_positions(0))*spline%slope_L
232  !lift(1,2) = -1._f64/3._f64*(spline%node_positions(nc)-spline%node_positions(nc-1))*spline%slope_R
233  !lift(2,1) = -2._f64/3._f64*(spline%node_positions(1)+spline%node_positions(2)-2.0_f64*spline%node_positions(0))*spline%slope_L
234  !lift(2,2) = 2._f64/3._f64*(2.0_f64*spline%node_positions(nc)-spline%node_positions(nc-1)-spline%node_positions(nc-2))*spline%slope_R
235 
236  lift(1, 1) = 1._f64/3._f64*(spline%node_positions(1) - spline%node_positions(0))*spline%slope_L
237  lift(1, 2) = -1._f64/3._f64*(spline%node_positions(nc) - spline%node_positions(nc - 1))*spline%slope_R
238  lift(2, 1) = -1._f64/3._f64*(spline%node_positions(1) - spline%node_positions(-1))&
239  &*(spline%node_positions(1) - spline%node_positions(-2))/(spline%node_positions(1) - spline%node_positions(0))*spline%slope_L
240  lift(2, 2) = 1._f64/3._f64*(spline%node_positions(nc + 1) - spline%node_positions(nc - 1))&
241  &*(spline%node_positions(nc+2)-spline%node_positions(nc-1))/(spline%node_positions(nc)-spline%node_positions(nc-1))*spline%slope_R
242 
243  lift(1:2, 1:2) = lift(1:2, 1:2)*(spline%xmax - spline%xmin)
244 
245  lift(3, 1) = (spline%node_positions(2) - spline%node_positions(0))*(spline%node_positions(1) - spline%node_positions(0))
246  lift(3,1) = lift(3,1)-(spline%node_positions(0)-spline%node_positions(-1))*(spline%node_positions(0)-spline%node_positions(-2))
247  lift(3,1) = lift(3,1)/((spline%node_positions(2)-spline%node_positions(-1))*(spline%node_positions(1)-spline%node_positions(0)))
248 
249  lift(3,2) = (spline%node_positions(nc+1)-spline%node_positions(nc))*(spline%node_positions(nc+2)-spline%node_positions(nc-1))
250  lift(3,2) = lift(3,2)/((spline%node_positions(nc+1)-spline%node_positions(nc-2))*(spline%node_positions(nc)-spline%node_positions(nc-1)))
251 
252  lift(4, 1) = (spline%node_positions(0) - spline%node_positions(-1))*(spline%node_positions(1) - spline%node_positions(-2))
253  lift(4,1) = lift(4,1)/((spline%node_positions(2)-spline%node_positions(-1))*(spline%node_positions(1)-spline%node_positions(0)))
254 
255 lift(4, 2) = (spline%node_positions(nc) - spline%node_positions(nc - 2))*(spline%node_positions(nc) - spline%node_positions(nc - 1))
256  lift(4,2) = lift(4,2)-(spline%node_positions(nc+2)-spline%node_positions(nc))*(spline%node_positions(nc+1)-spline%node_positions(nc))
257  lift(4,2) = lift(4,2)/((spline%node_positions(nc+1)-spline%node_positions(nc-2))*(spline%node_positions(nc)-spline%node_positions(nc-1)))
258 
259  call compute_spline_nonunif_1d_hermite_aux(fp, nc, buf, ibuf, coeffs, lift)
260  end subroutine compute_spline_nonunif_1d_hermite
261 
262 !0.95165885066842626 -0.66689806247771388
263 ! 0.16666666666666666 0.59333333333333371 0.23999999999999971
264 ! 3.07692307692324215E-003 0.20571428571428574 0.79120879120879106
265 ! 3.57142857142822359E-002 0.79761904761903923 0.16666666666667851
266 !
267 ! 0.16666666666667851 0.79761904761903923 3.57142857142822359E-002
268 ! 0.79120879120879106 0.20571428571428568 3.07692307692324129E-003
269 ! 0.23999999999999971 0.59333333333333360 0.16666666666666666
270 
272  subroutine sll_s_setup_spline_nonunif_1d_periodic_aux(node_pos, N, buf, ibuf)
273  sll_real64, dimension(:), pointer :: node_pos, buf
274  sll_int32, intent(in) :: n
275  sll_real64, dimension(:), pointer :: a, cts
276  sll_int32, dimension(:), pointer :: ipiv, ibuf
277  sll_int32 :: i
278  a => buf(1:3*n)
279  cts => buf(3*n + 1:10*n)
280  ipiv => ibuf(1:n)
281 
282  node_pos(-1) = node_pos(n - 1) - (node_pos(n) - node_pos(0))
283  node_pos(-2) = node_pos(n - 2) - (node_pos(n) - node_pos(0))
284  node_pos(n + 1) = node_pos(1) + (node_pos(n) - node_pos(0))
285  node_pos(n + 2) = node_pos(2) + (node_pos(n) - node_pos(0))
286 
287  !fill a with mesh information
288  do i = 0, n - 1
289  !subdiagonal terms
290  a(3*i+1)=(node_pos(i+1)-node_pos(i))*(node_pos(i+1)-node_pos(i))/((node_pos(i+1)-node_pos(i-1))*(node_pos(i+1)-node_pos(i-2)))
291  !a(3*i+1)=(node_pos(i+1)-node_pos(i))/(node_pos(i+1)-node_pos(i-1))*((node_pos(i+1)-node_pos(i))/(node_pos(i+1)-node_pos(i-2)))
292  !overdiagonal terms
293  !a(3*i+3)=(node_pos(i)-node_pos(i-1))*(node_pos(i)-node_pos(i-1))/((node_pos(i+1)-node_pos(i-1))*(node_pos(i+2)-node_pos(i-1)))
294  a(3*i+3)=(node_pos(i)-node_pos(i-1))/(node_pos(i+1)-node_pos(i-1))*((node_pos(i)-node_pos(i-1))/(node_pos(i+2)-node_pos(i-1)))
295  !diagonal terms
296  a(3*i + 2) = 1.0_f64 - a(3*i + 1) - a(3*i + 3)
297  end do
298  !print *,'#a=',a
299  !initialize the tridiagonal solver
300  call sll_s_setup_cyclic_tridiag(a, n, cts, ipiv)
301  !print *,'#cts=',maxval(cts),minval(cts)
302  !print *,'#ipiv=',ipiv
303  end subroutine sll_s_setup_spline_nonunif_1d_periodic_aux
304 
305  subroutine setup_spline_nonunif_1d_hermite_aux(node_pos, N, buf, ibuf)
306  sll_real64, dimension(:), pointer :: node_pos, buf
307  sll_int32, intent(in) :: n
308  sll_real64, dimension(:), pointer :: a, cts
309  sll_int32, dimension(:), pointer :: ipiv, ibuf
310  sll_int32 :: i, np
311  np = n + 1
312  a => buf(1:3*np)
313  cts => buf(3*np + 1:10*np)
314  ipiv => ibuf(1:np)
315 
316  !symmetric case
317  node_pos(-1) = 2._f64*node_pos(0) - node_pos(1)
318  node_pos(-2) = 2._f64*node_pos(0) - node_pos(2)
319  node_pos(n + 1) = 2.0_f64*node_pos(n) - node_pos(n - 1)
320  node_pos(n + 2) = 2.0_f64*node_pos(n) - node_pos(n - 2)
321 
322  !triple point case
323  node_pos(-1) = node_pos(0)
324  node_pos(-2) = node_pos(0)
325  node_pos(n + 1) = node_pos(n)
326  node_pos(n + 2) = node_pos(n)
327 
328  !fill a with mesh information
329  !a(1)=0.0_f64
330  !a(2)=(node_pos(2)-node_pos(0))/(node_pos(1)+node_pos(2)-2._f64*node_pos(0))
331  !a(3)=1.0_f64-a(2)
332  a(1) = 0.0_f64
333  a(2) = (node_pos(2) - node_pos(0))/(node_pos(2) - node_pos(-1))
334  a(3) = 1.0_f64 - a(2)
335  do i = 1, n - 1
336  !subdiagonal terms
337  a(3*i + 1) = (node_pos(i + 1) - node_pos(i))**2/((node_pos(i + 1) - node_pos(i - 1))*(node_pos(i + 1) - node_pos(i - 2)))
338  !overdiagonal terms
339  a(3*i + 3) = (node_pos(i) - node_pos(i - 1))**2/((node_pos(i + 1) - node_pos(i - 1))*(node_pos(i + 2) - node_pos(i - 1)))
340  !diagonal terms
341  a(3*i + 2) = 1.0_f64 - a(3*i + 1) - a(3*i + 3)
342  end do
343  !a(3*N+2)=(node_pos(N)-node_pos(N-2))/(2._f64*node_pos(N)-node_pos(N-1)-node_pos(N-2))
344  !a(3*N+1)=1.0_f64-a(3*N+2)
345  !a(3*N+3)=0.0_f64
346  a(3*n + 2) = (node_pos(n) - node_pos(n - 2))/(node_pos(n + 1) - node_pos(n - 2))
347  a(3*n + 1) = 1.0_f64 - a(3*n + 2)
348  a(3*n + 3) = 0.0_f64
349 
350  !initialize the tridiagonal solver
351  call sll_s_setup_cyclic_tridiag(a, np, cts, ipiv)
352 
353  end subroutine setup_spline_nonunif_1d_hermite_aux
354 
356  subroutine compute_spline_nonunif_1d_periodic_aux(f, N, buf, ibuf, coeffs)
357  sll_real64, dimension(:), pointer :: f, buf, coeffs
358  sll_int32, intent(in) :: n
359  sll_real64, dimension(:), pointer :: cts!, a
360  sll_int32, dimension(:), pointer :: ipiv, ibuf
361  !sll_real64 :: linf_err,tmp
362  !a => buf(1:3*N)
363  cts => buf(3*n + 1:10*n)
364  ipiv => ibuf(1:n)
365 
366  !compute the spline coefficients
367  call sll_o_solve_cyclic_tridiag(cts, ipiv, f, n, coeffs(0:n - 1))
368  coeffs(-1) = coeffs(n - 1)
369  coeffs(n) = coeffs(0)
370  coeffs(n + 1) = coeffs(1)
371 
372  !linf_err=0._f64
373  !do i=1,N
374  ! tmp=a(3*(i-1)+1)*coeffs(i-2)+a(3*(i-1)+2)*coeffs(i-1)+a(3*(i-1)+3)*coeffs(i)-f(i)
375  ! if(abs(tmp)>linf_err)then
376  ! linf_err=abs(tmp)
377  ! endif
378  !enddo
379  !print *,'error of compute_spline=',linf_err
380  end subroutine compute_spline_nonunif_1d_periodic_aux
381 
383  subroutine sll_s_compute_spline_nonunif_1d_periodic_aux2(f, N, buf, ibuf, coeffs)
384  sll_real64, dimension(:) :: f
385  sll_real64, dimension(:), pointer :: buf
386  sll_real64, dimension(:), pointer :: coeffs
387  sll_int32, intent(in) :: n
388  sll_real64, dimension(:), pointer :: cts!, a
389  sll_int32, dimension(:), pointer :: ipiv, ibuf
390  sll_real64, dimension(:), pointer :: a
391  sll_real64, dimension(:), allocatable :: b
392  sll_real64, dimension(:), pointer :: x
393  sll_int32 :: ierr
394 
395  !sll_real64 :: linf_err,tmp
396  !a => buf(1:3*N)
397  cts => buf(3*n + 1:10*n)
398  ipiv => ibuf(1:n)
399 
400  a => buf(1:3*n)
401  !b => f(1:N)
402  x => coeffs(0:n - 1)
403 
404  sll_allocate(b(n), ierr)
405  b(1:n) = f(1:n)
406 
407  !compute the spline coefficients
408  !call sll_o_solve_cyclic_tridiag( cts, ipiv, f, N, coeffs(0:N-1) )
409  coeffs(0:n - 1) = 0._f64
410  !do i=0,N-1
411  ! if(coeffs(i) .ne. 0._f64)then
412  ! print *,'#coeffs=',coeffs
413  ! stop
414  ! endif
415  !enddo
416 
417  call sll_o_solve_cyclic_tridiag(cts, ipiv, f(1:n), n, coeffs(0:n - 1))
418 
419  coeffs(-1) = coeffs(n - 1)
420  coeffs(n) = coeffs(0)
421  coeffs(n + 1) = coeffs(1)
422 
423  !linf_err=0._f64
424  !do i=1,N
425  ! tmp=a(3*(i-1)+1)*coeffs(i-2)+a(3*(i-1)+2)*coeffs(i-1)+a(3*(i-1)+3)*coeffs(i)-f(i)
426  ! if(abs(tmp)>linf_err)then
427  ! linf_err=abs(tmp)
428  ! endif
429  !enddo
430  !print *,'error of compute_spline=',linf_err
431  end subroutine sll_s_compute_spline_nonunif_1d_periodic_aux2
432 
433  subroutine compute_spline_nonunif_1d_hermite_aux(f, N, buf, ibuf, coeffs, lift)
434  sll_real64, dimension(:), pointer :: f, buf, coeffs
435  sll_int32, intent(in) :: n
436  sll_real64, intent(in) :: lift(4, 2)
437  sll_real64, dimension(:), pointer :: cts
438  sll_int32, dimension(:), pointer :: ipiv, ibuf
439  sll_int32 :: np
440  np = n + 1
441  cts => buf(3*np + 1:10*np)
442  ipiv => ibuf(1:np)
443 
444  !compute the spline coefficients with inplace tridiagonal solver
445  coeffs(1:n - 1) = f(2:n)
446  coeffs(0) = f(1) + lift(1, 1)
447  coeffs(n) = f(n + 1) + lift(1, 2)
448  call sll_o_solve_cyclic_tridiag(cts, ipiv, coeffs(0:n), np, coeffs(0:n))
449  !coeffs(-1) = coeffs(1)+lift(2,1)
450  !coeffs(N+1) = coeffs(N-1)+lift(2,2)
451 
452  coeffs(-1) = lift(3, 1)*coeffs(0) + lift(4, 1)*coeffs(1) + lift(2, 1)
453  coeffs(n + 1) = lift(3, 2)*coeffs(n - 1) + lift(4, 2)*coeffs(n) + lift(2, 2)
454 
455  end subroutine compute_spline_nonunif_1d_hermite_aux
456 
458  function interpolate_value_nonunif(x, spline)
459  !warning not tested, should not work->t o be adapted like array
460  sll_real64 :: interpolate_value_nonunif
461  sll_real64, intent(in) :: x
462  type(sll_t_cubic_nonunif_spline_1d), pointer :: spline
463  sll_int32 :: j
464  sll_real64 ::xx
465  sll_real64, dimension(:), pointer :: xj, coef
466  sll_real64 :: w(0:3)
467  xx = (x - spline%xmin)/(spline%xmax - spline%xmin)
468  if (.not. ((xx .ge. 0.0_f64) .and. (xx .le. 1.0_f64))) then
469  print *, 'bad_value of x=', x, 'xmin=', spline%xmin, 'xmax=', spline%xmax, xx, xx - 1.0_f64
470  print *, 'in function interpolate_value_nonunif()'
471  stop
472  end if
473  !SLL_ASSERT( (xx .ge. 0.0_f64) .and. (xx .le. 1.0_f64 ))
474 
475  !localization of xx
476  if (xx == 1.0_f64) then
477  j = spline%n_cells
478  else
479  j = 0
480  do while (spline%node_positions(j) .le. xx)
481  j = j + 1
482  end do
483  end if
484  xj => spline%node_positions(j:)
485 
486  sll_assert((xx .ge. xj(0)) .and. (xx .le. xj(1)))
487  sll_assert((xj(0) == spline%node_positions(j - 1)) .and. (xj(1) == spline%node_positions(j)))
488 
489  !compute weights
490  w(0) = (xj(1) - xx)*(xj(1) - xx)*(xj(1) - xx)/((xj(1) - xj(0))*(xj(1) - xj(-1))*(xj(1) - xj(-2)));
491  w(1) = (xj(1) - xx)*(xj(1) - xx)*(xx - xj(-2))/((xj(1) - xj(0))*(xj(1) - xj(-1))*(xj(1) - xj(-2)));
492  w(1) = w(1) + (xj(2) - xx)*(xj(1) - xx)*(xx - xj(-1))/((xj(1) - xj(0))*(xj(1) - xj(-1))*(xj(2) - xj(-1)));
493  w(1) = w(1) + (xj(2) - xx)*(xj(2) - xx)*(xx - xj(0))/((xj(1) - xj(0))*(xj(2) - xj(0))*(xj(2) - xj(-1)));
494  w(2) = (xj(1) - xx)*(xx - xj(-1))*(xx - xj(-1))/((xj(1) - xj(0))*(xj(1) - xj(-1))*(xj(2) - xj(-1)));
495  w(2) = w(2) + (xj(2) - xx)*(xx - xj(-1))*(xx - xj(0))/((xj(1) - xj(0))*(xj(2) - xj(0))*(xj(2) - xj(-1)));
496  w(2) = w(2) + (xj(3) - xx)*(xx - xj(0))*(xx - xj(0))/((xj(1) - xj(0))*(xj(2) - xj(0))*(xj(3) - xj(0)));
497  w(3) = (xx - xj(0))*(xx - xj(0))*(xx - xj(0))/((xj(1) - xj(0))*(xj(2) - xj(0))*(xj(3) - xj(0)));
498  coef => spline%coeffs(j - 1:)
499 
500  interpolate_value_nonunif = w(0)*coef(0) + w(1)*coef(1) + w(2)*coef(2) + w(3)*coef(3)
501  end function interpolate_value_nonunif
502 
504  subroutine sll_s_interpolate_array_value_nonunif(a_in, a_out, n, spline)
505  sll_int32, intent(in) :: n
506  sll_real64, dimension(1:n), intent(in) :: a_in
507  sll_real64, dimension(1:n), intent(out) :: a_out
508  sll_real64 :: x
509  type(sll_t_cubic_nonunif_spline_1d), pointer :: spline
510  sll_int32 :: i, j, shift = 3
511  sll_real64 ::xx
512  sll_real64, dimension(:), pointer :: xj
513  sll_real64, dimension(:), pointer :: coef
514  sll_real64 :: w(4)
515 
516  x = a_in(1)
517  xx = (x - spline%xmin)/(spline%xmax - spline%xmin)
518  if (.not. ((xx .ge. 0.0_f64) .and. (xx .le. 1.0_f64))) then
519  print *, 'bad_value of x=', x, 'xmin=', spline%xmin, 'xmax=', spline%xmax, xx, xx - 1.0_f64
520  print *, 'in subroutine sll_s_interpolate_array_value_nonunif()'
521  stop
522  end if
523  !SLL_ASSERT( (xx .ge. 0.0_f64) .and. (xx .le. 1.0_f64 ))
524 
525  !localization of xx
526  j = 0
527  if (xx == 1.0_f64) then
528  j = spline%n_cells
529  else
530  do while (spline%node_positions(j) .le. xx)
531  j = j + 1
532  end do
533  end if
534 
535  do i = 1, n
536 
537  x = (a_in(i) - spline%xmin)/(spline%xmax - spline%xmin)
538  if (.not. ((x .ge. 0.0_f64) .and. (x .le. 1.0_f64))) then
539  print *, 'bad_value of a_in(', i, ')=', a_in(i), 'xmin=', spline%xmin, 'xmax=', spline%xmax
540  print *, 'in subroutine sll_s_interpolate_array_value_nonunif()'
541  stop
542  end if
543  if (x == 1.0_f64) then
544  j = spline%n_cells
545  else
546  if (x .ge. xx) then
547  do while (spline%node_positions(j) .le. x)
548  j = j + 1
549  end do
550  else
551  do while (spline%node_positions(j) .gt. x)
552  j = j - 1
553  end do
554  j = j + 1
555  end if
556  end if
557  xx = x
558  xj => spline%node_positions(j - shift:)
559 
560  !print *,i,Xj(0+shift),xx,Xj(1+shift)
561  !stop
562  if (.not. ((xx .ge. xj(0 + shift)) .and. (xx .lt. xj(1 + shift)))) then
563  if (xx .ne. 1.0_f64) then
564  print *, xj(0 + shift), xx, xj(1 + shift)
565 
566  stop
567  end if
568  end if
569  sll_assert((xx .ge. xj(0 + shift)) .and. (xx .le. xj(1 + shift)))
570  sll_assert((xj(0 + shift) == spline%node_positions(j - 1)) .and. (xj(1 + shift) == spline%node_positions(j)))
571 
572  !compute weights
573  w(1)=(xj(shift+1)-xx)*(xj(shift+1)-xx)*(xj(shift+1)-xx)/((xj(shift+1)-xj(shift+0))*(xj(shift+1)-xj(shift-1))*(xj(shift+1)-xj(shift-2)));
574  w(2)=(xj(shift+1)-xx)*(xj(shift+1)-xx)*(xx-xj(shift-2))/((xj(shift+1)-xj(shift+0))*(xj(shift+1)-xj(shift-1))*(xj(shift+1)-xj(shift-2)));
575  w(2)=w(2)+(xj(shift+2)-xx)*(xj(shift+1)-xx)*(xx-xj(shift-1))/((xj(shift+1)-xj(shift+0))*(xj(shift+1)-xj(shift-1))*(xj(shift+2)-xj(shift-1)));
576  w(2)=w(2)+(xj(shift+2)-xx)*(xj(shift+2)-xx)*(xx-xj(shift+0))/((xj(shift+1)-xj(shift+0))*(xj(shift+2)-xj(shift+0))*(xj(shift+2)-xj(shift-1)));
577  w(3)=(xj(shift+1)-xx)*(xx-xj(shift-1))*(xx-xj(shift-1))/((xj(shift+1)-xj(shift+0))*(xj(shift+1)-xj(shift-1))*(xj(shift+2)-xj(shift-1)));
578  w(3)=w(3)+(xj(shift+2)-xx)*(xx-xj(shift-1))*(xx-xj(shift+0))/((xj(shift+1)-xj(shift+0))*(xj(shift+2)-xj(shift+0))*(xj(shift+2)-xj(shift-1)));
579  w(3)=w(3)+(xj(shift+3)-xx)*(xx-xj(shift+0))*(xx-xj(shift+0))/((xj(shift+1)-xj(shift+0))*(xj(shift+2)-xj(shift+0))*(xj(shift+3)-xj(shift+0)));
580  w(4)=(xx-xj(shift+0))*(xx-xj(shift+0))*(xx-xj(shift+0))/((xj(shift+1)-xj(shift+0))*(xj(shift+2)-xj(shift+0))*(xj(shift+3)-xj(shift+0)));
581  !coef => spline%coeffs(j-1:)
582  !print *,i,xx,j,w(0),w(1),w(2),w(3),Xj(-2:3)
583  !a_out(i) = w(0)*coef(0)+w(1)*coef(1)+w(2)*coef(2)+w(3)*coef(3)
584 
585  coef => spline%coeffs(j - 2:)
586  !print *,i,xx,j,w(0),w(1),w(2),w(3),Xj(-2:3)
587  a_out(i) = w(1)*coef(1) + w(2)*coef(2) + w(3)*coef(3) + w(4)*coef(4)
588 
589  end do
590  end subroutine sll_s_interpolate_array_value_nonunif
591 
593  subroutine sll_s_interpolate_array_value_nonunif_aux(a_in, a_out, n, node_pos, coeffs, n_cells)
594  sll_int32, intent(in) :: n, n_cells
595  sll_real64, dimension(1:n), intent(in) :: a_in
596  !sll_real64, dimension(-1:n+1), intent(in) :: node_pos
597  sll_real64, dimension(1:n), intent(out) :: a_out
598  sll_real64 :: x
599  !type(sll_t_cubic_nonunif_spline_1d), pointer :: spline
600  sll_int32 :: i, j, shift = 3
601  sll_real64 ::xx
602  sll_real64, dimension(:), pointer :: xj
603  sll_real64, dimension(:), pointer :: coef, coeffs, node_pos
604  sll_real64 :: w(4)
605 
606  !do i=-1,n_cells+1
607  ! print *,i,node_pos(i)
608  !enddo
609 
610  xx = a_in(1)
611  !xx = (x-spline%xmin)/(spline%xmax-spline%xmin)
612  if (.not. ((xx .ge. 0.0_f64) .and. (xx .le. 1.0_f64))) then
613  print *, 'bad_value of x=', x!, 'xmin=', spline%xmin, 'xmax=', spline%xmax, xx,xx-1.0_f64
614  print *, 'in subroutine sll_s_interpolate_array_value_nonunif()'
615  stop
616  end if
617  !SLL_ASSERT( (xx .ge. 0.0_f64) .and. (xx .le. 1.0_f64 ))
618 
619  !localization of xx
620  j = 0
621  if (xx == 1.0_f64) then
622  j = n_cells
623  else
624  do while (node_pos(j) .le. xx)
625  j = j + 1
626  end do
627  end if
628 
629  do i = 1, n
630 
631  x = a_in(i)
632  if (.not. ((x .ge. 0.0_f64) .and. (x .le. 1.0_f64))) then
633  print *, 'bad_value of a_in(', i, ')=', a_in(i)!, 'xmin=', spline%xmin, 'xmax=', spline%xmax
634  print *, 'in subroutine sll_s_interpolate_array_value_nonunif()'
635  stop
636  end if
637  if (x == 1.0_f64) then
638  j = n_cells!spline%n_cells
639  else
640  if (x .ge. xx) then
641  do while (node_pos(j) .le. x)
642  j = j + 1
643  end do
644  else
645  do while (node_pos(j) .gt. x)
646  j = j - 1
647  end do
648  j = j + 1
649  end if
650  end if
651  xx = x
652  xj => node_pos(j - shift:)
653 
654  !print *,i,Xj(0+shift),xx,Xj(1+shift)
655  !print *,i,Xj(-2+shift:3+shift)
656  !stop
657  if (.not. ((xx .ge. xj(0 + shift)) .and. (xx .lt. xj(1 + shift)))) then
658  if (xx .ne. 1.0_f64) then
659  print *, xj(0 + shift), xx, xj(1 + shift)
660 
661  stop
662  end if
663  end if
664  !SLL_ASSERT( (xx .ge. Xj(0+shift)) .and. (xx .le. Xj(1+shift) ))
665  !SLL_ASSERT( (Xj(0+shift)==spline%node_positions(j-1)) .and. (Xj(1+shift)==spline%node_positions(j)))
666 
667  !compute weights
668  w(1) = (xj(shift + 1) - xx)*(xj(shift + 1) - xx)*(xj(shift + 1) - xx)&
669  &/((xj(shift + 1) - xj(shift + 0))*(xj(shift + 1) - xj(shift - 1))*(xj(shift + 1) - xj(shift - 2)));
670  w(2) = (xj(shift + 1) - xx)*(xj(shift + 1) - xx)*(xx - xj(shift - 2))&
671  &/((xj(shift + 1) - xj(shift + 0))*(xj(shift + 1) - xj(shift - 1))*(xj(shift + 1) - xj(shift - 2)));
672  w(2) = w(2) + (xj(shift + 2) - xx)*(xj(shift + 1) - xx)*(xx - xj(shift - 1))&
673  &/((xj(shift + 1) - xj(shift + 0))*(xj(shift + 1) - xj(shift - 1))*(xj(shift + 2) - xj(shift - 1)));
674  w(2) = w(2) + (xj(shift + 2) - xx)*(xj(shift + 2) - xx)*(xx - xj(shift + 0))&
675  &/((xj(shift + 1) - xj(shift + 0))*(xj(shift + 2) - xj(shift + 0))*(xj(shift + 2) - xj(shift - 1)));
676  w(3) = (xj(shift + 1) - xx)*(xx - xj(shift - 1))*(xx - xj(shift - 1))&
677  &/((xj(shift + 1) - xj(shift + 0))*(xj(shift + 1) - xj(shift - 1))*(xj(shift + 2) - xj(shift - 1)));
678  w(3) = w(3) + (xj(shift + 2) - xx)*(xx - xj(shift - 1))*(xx - xj(shift + 0))&
679  &/((xj(shift + 1) - xj(shift + 0))*(xj(shift + 2) - xj(shift + 0))*(xj(shift + 2) - xj(shift - 1)));
680  w(3) = w(3) + (xj(shift + 3) - xx)*(xx - xj(shift + 0))*(xx - xj(shift + 0))&
681  &/((xj(shift + 1) - xj(shift + 0))*(xj(shift + 2) - xj(shift + 0))*(xj(shift + 3) - xj(shift + 0)));
682  w(4) = (xx - xj(shift + 0))*(xx - xj(shift + 0))*(xx - xj(shift + 0))&
683  &/((xj(shift + 1) - xj(shift + 0))*(xj(shift + 2) - xj(shift + 0))*(xj(shift + 3) - xj(shift + 0)));
684  !coef => spline%coeffs(j-1:)
685  !print *,i,xx,j,w(1),w(2),w(3),w(4)!,Xj(-2:3)
686  !a_out(i) = w(0)*coef(0)+w(1)*coef(1)+w(2)*coef(2)+w(3)*coef(3)
687 
688  coef => coeffs(j - 2:)
689  a_out(i) = w(1)*coef(1) + w(2)*coef(2) + w(3)*coef(3) + w(4)*coef(4)
690  !print *,i,xx,j,w(1),w(2),w(3),w(4),coef(1:4),a_out(i)
691  !print *,Xj(shift-2:shift+3)
692  end do
693  !do i=1,n
694  ! print *,i,a_in(i),a_out(i)
695  !enddo
696  end subroutine sll_s_interpolate_array_value_nonunif_aux
697 
698 end module sll_m_cubic_non_uniform_splines
699 
sll_m_cubic_non_uniform_splines::sll_o_delete
Definition: sll_m_cubic_non_uniform_splines.F90:57
sll_m_tridiagonal::sll_o_solve_cyclic_tridiag
Solve tridiagonal system (double or complex)
Definition: sll_m_tridiagonal.F90:72
sll_m_boundary_condition_descriptors
Describe different boundary conditions.
Definition: sll_m_boundary_condition_descriptors.F90:20
sll_m_cubic_non_uniform_splines
Provides capabilities for data interpolation with cubic B-splines on non-uniform meshes.
Definition: sll_m_cubic_non_uniform_splines.F90:6
sll_m_cubic_non_uniform_splines::sll_s_compute_spline_nonunif
subroutine, public sll_s_compute_spline_nonunif(f, spline, node_positions, sl, sr)
compute splines coefficients
Definition: sll_m_cubic_non_uniform_splines.F90:115
sll_m_cubic_non_uniform_splines::sll_s_cubic_nonunif_spline_1d_init
subroutine sll_s_cubic_nonunif_spline_1d_init(self, n_cells, bc_type)
Initialize spline object.
Definition: sll_m_cubic_non_uniform_splines.F90:78
sll_m_cubic_non_uniform_splines::sll_f_new_cubic_nonunif_spline_1d
type(sll_t_cubic_nonunif_spline_1d) function, pointer, public sll_f_new_cubic_nonunif_spline_1d(n_cells, bc_type)
create new spline object
Definition: sll_m_cubic_non_uniform_splines.F90:65
sll_m_cubic_non_uniform_splines::setup_spline_nonunif_1d_hermite_aux
subroutine setup_spline_nonunif_1d_hermite_aux(node_pos, N, buf, ibuf)
Definition: sll_m_cubic_non_uniform_splines.F90:306
sll_m_cubic_non_uniform_splines::interpolate_value_nonunif
real(kind=f64) function interpolate_value_nonunif(x, spline)
get spline interpolate at point x
Definition: sll_m_cubic_non_uniform_splines.F90:459
sll_m_cubic_non_uniform_splines::sll_s_interpolate_array_value_nonunif_aux
subroutine, public sll_s_interpolate_array_value_nonunif_aux(a_in, a_out, n, node_pos, coeffs, n_cells)
Definition: sll_m_cubic_non_uniform_splines.F90:594
sll_m_cubic_non_uniform_splines::delete_cubic_nonunif_spline_1d
subroutine delete_cubic_nonunif_spline_1d(spline, ierr)
delete spline object
Definition: sll_m_cubic_non_uniform_splines.F90:103
sll_m_cubic_non_uniform_splines::compute_spline_nonunif_1d_periodic_aux
subroutine compute_spline_nonunif_1d_periodic_aux(f, N, buf, ibuf, coeffs)
Definition: sll_m_cubic_non_uniform_splines.F90:357
sll_m_cubic_non_uniform_splines::sll_s_interpolate_array_value_nonunif
subroutine, public sll_s_interpolate_array_value_nonunif(a_in, a_out, n, spline)
get spline interpolate on an array of points
Definition: sll_m_cubic_non_uniform_splines.F90:505
sll_m_cubic_non_uniform_splines::compute_spline_nonunif_1d_periodic
subroutine compute_spline_nonunif_1d_periodic(f, spline)
Definition: sll_m_cubic_non_uniform_splines.F90:176
sll_m_cubic_non_uniform_splines::compute_spline_nonunif_1d_hermite
subroutine compute_spline_nonunif_1d_hermite(f, spline)
Definition: sll_m_cubic_non_uniform_splines.F90:205
sll_m_cubic_non_uniform_splines::sll_s_compute_spline_nonunif_1d_periodic_aux2
subroutine, public sll_s_compute_spline_nonunif_1d_periodic_aux2(f, N, buf, ibuf, coeffs)
Definition: sll_m_cubic_non_uniform_splines.F90:384
sll_m_cubic_non_uniform_splines::sll_s_setup_spline_nonunif_1d_periodic_aux
subroutine, public sll_s_setup_spline_nonunif_1d_periodic_aux(node_pos, N, buf, ibuf)
Definition: sll_m_cubic_non_uniform_splines.F90:273
sll_m_cubic_non_uniform_splines::compute_spline_nonunif_1d_hermite_aux
subroutine compute_spline_nonunif_1d_hermite_aux(f, N, buf, ibuf, coeffs, lift)
Definition: sll_m_cubic_non_uniform_splines.F90:434
sll_m_tridiagonal
Tridiagonal system solver.
Definition: sll_m_tridiagonal.F90:57
sll_m_tridiagonal::sll_s_setup_cyclic_tridiag
subroutine, public sll_s_setup_cyclic_tridiag(a, n, cts, ipiv)
Give the factorization of the matrix in argument.
Definition: sll_m_tridiagonal.F90:162
sll_m_cubic_non_uniform_splines::sll_t_cubic_nonunif_spline_1d
Spline object.
Definition: sll_m_cubic_non_uniform_splines.F90:36
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:41 for Semi-Lagrangian Library by   doxygen 1.9.1