Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_cubic_spline_interpolator_2d.F90
Go to the documentation of this file.
1 !**************************************************************
2 ! Copyright INRIA
3 ! Authors :
4 ! CALVI project team
5 !
6 ! This code SeLaLib (for Semi-Lagrangian-Library)
7 ! is a parallel library for simulating the plasma turbulence
8 ! in a tokamak.
9 !
10 ! This software is governed by the CeCILL-B license
11 ! under French law and abiding by the rules of distribution
12 ! of free software. You can use, modify and redistribute
13 ! the software under the terms of the CeCILL-B license as
14 ! circulated by CEA, CNRS and INRIA at the following URL
15 ! "http://www.cecill.info".
16 !**************************************************************
17 
24 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
25 #include "sll_assert.h"
26 #include "sll_errors.h"
27 #include "sll_memory.h"
28 #include "sll_working_precision.h"
29 
31  sll_p_hermite, &
32  sll_p_periodic
33 
34  use sll_m_cubic_splines, only: &
36  sll_f_cubic_spline_2d_get_x1_delta, &
37  sll_f_cubic_spline_2d_get_x1_max, &
38  sll_f_cubic_spline_2d_get_x1_min, &
39  sll_f_cubic_spline_2d_get_x2_delta, &
40  sll_f_cubic_spline_2d_get_x2_max, &
41  sll_f_cubic_spline_2d_get_x2_min, &
48 
49  use sll_m_interpolators_2d_base, only: &
51 
52  implicit none
53 
54  public :: &
58 
59  private
60 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
61 
71  sll_int32 :: npts1
73  sll_int32 :: npts2
75  type(sll_t_cubic_spline_2d) :: spline
77  sll_int32 :: bc_type1
79  sll_int32 :: bc_type2
81  sll_real64, dimension(:, :), pointer :: interpolation_points
82  contains
84  procedure, pass(interpolator) :: init => initialize_cs2d_interpolator
86  procedure :: compute_interpolants => compute_interpolants_cs2d
88  procedure :: interpolate_from_interpolant_value => interpolate_value_cs2d
90  procedure :: interpolate_from_interpolant_derivative_eta1 => interpolate_deriv1_cs2d
92  procedure :: interpolate_from_interpolant_derivative_eta2 => interpolate_deriv2_cs2d
94  procedure, pass :: interpolate_array => spline_interpolate2d
96  procedure, pass :: interpolate_array_disp => spline_interpolate2d_disp
98  procedure, pass :: set_coefficients => set_coefficients_cs2d
100  procedure, pass :: get_coefficients => get_coefficients_cs2d
102  procedure, pass :: coefficients_are_set => coefficients_are_set_cs2d
104  procedure, pass :: delete => delete_sll_cubic_spline_interpolator_2d
106 
109  type(sll_t_cubic_spline_interpolator_2d), pointer :: interp
111 
113  interface sll_o_delete
115  end interface sll_o_delete
116 
117 contains
118 
120  class(sll_t_cubic_spline_interpolator_2d), intent(inout) :: interpolator
121  call sll_s_cubic_spline_2d_free(interpolator%spline)
123 
128  npts1, &
129  npts2, &
130  eta1_min, &
131  eta1_max, &
132  eta2_min, &
133  eta2_max, &
134  eta1_bc_type, &
135  eta2_bc_type, &
136  const_eta1_min_slope, &
137  const_eta1_max_slope, &
138  const_eta2_min_slope, &
139  const_eta2_max_slope, &
140  eta1_min_slopes, &
141  eta1_max_slopes, &
142  eta2_min_slopes, &
143  eta2_max_slopes) &
144  result(interpolator)
145 
146  type(sll_t_cubic_spline_interpolator_2d), pointer :: interpolator
147 
148  sll_int32, intent(in) :: npts1
149  sll_int32, intent(in) :: npts2
150  sll_real64, intent(in) :: eta1_min
151  sll_real64, intent(in) :: eta1_max
152  sll_real64, intent(in) :: eta2_min
153  sll_real64, intent(in) :: eta2_max
154  sll_int32, intent(in) :: eta1_bc_type
155  sll_int32, intent(in) :: eta2_bc_type
156  sll_real64, intent(in), optional :: const_eta1_min_slope
157  sll_real64, intent(in), optional :: const_eta1_max_slope
158  sll_real64, intent(in), optional :: const_eta2_min_slope
159  sll_real64, intent(in), optional :: const_eta2_max_slope
160  sll_real64, dimension(:), intent(in), optional :: eta1_min_slopes
161  sll_real64, dimension(:), intent(in), optional :: eta1_max_slopes
162  sll_real64, dimension(:), intent(in), optional :: eta2_min_slopes
163  sll_real64, dimension(:), intent(in), optional :: eta2_max_slopes
164  sll_int32 :: ierr
165 
166  sll_allocate(interpolator, ierr)
167 
168  call interpolator%init( &
169  npts1, &
170  npts2, &
171  eta1_min, &
172  eta1_max, &
173  eta2_min, &
174  eta2_max, &
175  eta1_bc_type, &
176  eta2_bc_type, &
177  const_eta1_min_slope, &
178  const_eta1_max_slope, &
179  const_eta2_min_slope, &
180  const_eta2_max_slope, &
181  eta1_min_slopes, &
182  eta1_max_slopes, &
183  eta2_min_slopes, &
184  eta2_max_slopes)
185 
187 
188  ! We allow to use the enumerators of the splines module in this interpolator
189  ! because:
190  ! a. This is just a wrapper and is intimately related to the underlying
191  ! cubic splines module.
192  ! b. There is no uniform interface for the initialization anyway.
193  ! The underlying implementation with the splines module could be hidden but
194  ! I can't see a compelling reason why.
196  interpolator, &
197  npts1, &
198  npts2, &
199  eta1_min, &
200  eta1_max, &
201  eta2_min, &
202  eta2_max, &
203  eta1_bc_type, &
204  eta2_bc_type, &
205  const_eta1_min_slope, &
206  const_eta1_max_slope, &
207  const_eta2_min_slope, &
208  const_eta2_max_slope, &
209  eta1_min_slopes, &
210  eta1_max_slopes, &
211  eta2_min_slopes, &
212  eta2_max_slopes)
213 
214  class(sll_t_cubic_spline_interpolator_2d), intent(inout) :: interpolator
215  sll_int32, intent(in) :: npts1
216  sll_int32, intent(in) :: npts2
217  sll_real64, intent(in) :: eta1_min
218  sll_real64, intent(in) :: eta1_max
219  sll_real64, intent(in) :: eta2_min
220  sll_real64, intent(in) :: eta2_max
221  sll_int32, intent(in) :: eta1_bc_type
222  sll_int32, intent(in) :: eta2_bc_type
223  sll_real64, intent(in), optional :: const_eta1_min_slope
224  sll_real64, intent(in), optional :: const_eta1_max_slope
225  sll_real64, intent(in), optional :: const_eta2_min_slope
226  sll_real64, intent(in), optional :: const_eta2_max_slope
227  sll_real64, dimension(:), intent(in), optional :: eta1_min_slopes
228  sll_real64, dimension(:), intent(in), optional :: eta1_max_slopes
229  sll_real64, dimension(:), intent(in), optional :: eta2_min_slopes
230  sll_real64, dimension(:), intent(in), optional :: eta2_max_slopes
231 
232  interpolator%npts1 = npts1
233  interpolator%npts2 = npts2
234  interpolator%bc_type1 = eta1_bc_type
235  interpolator%bc_type2 = eta2_bc_type
236 
238  interpolator%spline, &
239  npts1, &
240  npts2, &
241  eta1_min, &
242  eta1_max, &
243  eta2_min, &
244  eta2_max, &
245  eta1_bc_type, &
246  eta2_bc_type, &
247  const_slope_x1_min=const_eta1_min_slope, &
248  const_slope_x1_max=const_eta1_max_slope, &
249  const_slope_x2_min=const_eta2_min_slope, &
250  const_slope_x2_max=const_eta2_max_slope, &
251  x1_min_slopes=eta1_min_slopes, &
252  x1_max_slopes=eta1_max_slopes, &
253  x2_min_slopes=eta2_min_slopes, &
254  x2_max_slopes=eta2_max_slopes)
255 
256  end subroutine
257 
259  interpolator, &
260  data_array, &
261  eta1_coords, &
262  size_eta1_coords, &
263  eta2_coords, &
264  size_eta2_coords)
265  class(sll_t_cubic_spline_interpolator_2d), intent(inout) :: interpolator
266  sll_real64, dimension(:, :), intent(in) :: data_array
267  sll_real64, dimension(:), intent(in), optional :: eta1_coords
268  sll_real64, dimension(:), intent(in), optional :: eta2_coords
269  sll_int32, intent(in), optional :: size_eta1_coords
270  sll_int32, intent(in), optional :: size_eta2_coords
271 
272  if (present(eta1_coords) .or. present(eta2_coords) &
273  .or. present(size_eta1_coords) .or. present(size_eta2_coords)) then
274  sll_error('compute_interpolants_cs2d', 'This case is not yet implemented')
275  end if
276 
277  call sll_s_cubic_spline_2d_compute_interpolant(data_array, interpolator%spline)
278 
279  end subroutine
280 
281  function interpolate_value_cs2d(interpolator, eta1, eta2) result(val)
282  class(sll_t_cubic_spline_interpolator_2d), intent(in) :: interpolator
283  sll_real64 :: val
284  sll_real64, intent(in) :: eta1
285  sll_real64, intent(in) :: eta2
286  val = sll_f_cubic_spline_2d_eval(eta1, eta2, interpolator%spline)
287  end function
288 
289  function interpolate_deriv1_cs2d(interpolator, eta1, eta2) result(val)
290  class(sll_t_cubic_spline_interpolator_2d), intent(in) :: interpolator
291  sll_real64 :: val
292  sll_real64, intent(in) :: eta1
293  sll_real64, intent(in) :: eta2
294  val = sll_f_cubic_spline_2d_eval_deriv_x1(eta1, eta2, interpolator%spline)
295  end function
296 
297  function interpolate_deriv2_cs2d(interpolator, eta1, eta2) result(val)
298  class(sll_t_cubic_spline_interpolator_2d), intent(in) :: interpolator
299  sll_real64 :: val
300  sll_real64, intent(in) :: eta1
301  sll_real64, intent(in) :: eta2
302 
303  val = sll_f_cubic_spline_2d_eval_deriv_x2(eta1, eta2, interpolator%spline)
304 
305  end function
306 
307  subroutine spline_interpolate2d(this, num_points1, num_points2, data_in, &
308  eta1, eta2, data_out)
309 
310  class(sll_t_cubic_spline_interpolator_2d), intent(inout) :: this
311  sll_int32, intent(in) :: num_points1
312  sll_int32, intent(in) :: num_points2
313  sll_real64, dimension(:, :), intent(in) :: eta1
314  sll_real64, dimension(:, :), intent(in) :: eta2
315  sll_real64, dimension(:, :), intent(in) :: data_in
316  sll_real64, intent(out) :: data_out(num_points1, num_points2)
317  ! local variables
318  sll_int32 :: i, j
319  ! compute the interpolating spline coefficients
320  call sll_s_cubic_spline_2d_compute_interpolant(data_in, this%spline)
321  do j = 1, num_points2
322  do i = 1, num_points1
323  data_out(i, j) = this%interpolate_from_interpolant_value(eta1(i, j), eta2(i, j))
324  end do
325  end do
326 
327  end subroutine spline_interpolate2d
328 
329  subroutine spline_interpolate2d_disp(this, &
330  num_points1, &
331  num_points2, &
332  data_in, &
333  alpha1, &
334  alpha2, &
335  data_out)
336 
337  class(sll_t_cubic_spline_interpolator_2d), intent(inout) :: this
338 
339  sll_int32, intent(in) :: num_points1
340  sll_int32, intent(in) :: num_points2
341  sll_real64, dimension(:, :), intent(in) :: alpha1
342  sll_real64, dimension(:, :), intent(in) :: alpha2
343  sll_real64, dimension(:, :), intent(in) :: data_in
344  sll_real64, intent(out) :: data_out(num_points1, num_points2)
345  sll_real64 :: eta1
346  sll_real64 :: eta1_min
347  sll_real64 :: eta1_max
348  sll_real64 :: delta_eta1
349  sll_real64 :: eta2
350  sll_real64 :: eta2_min
351  sll_real64 :: eta2_max
352  sll_real64 :: delta_eta2
353  sll_int32 :: i
354  sll_int32 :: j
355 
356  eta1_min = sll_f_cubic_spline_2d_get_x1_min(this%spline) !this%spline%x1_min
357  eta1_max = sll_f_cubic_spline_2d_get_x1_max(this%spline) !this%spline%x1_max
358  eta2_min = sll_f_cubic_spline_2d_get_x2_min(this%spline) !this%spline%x2_min
359  eta2_max = sll_f_cubic_spline_2d_get_x2_max(this%spline) !this%spline%x2_max
360  delta_eta1 = sll_f_cubic_spline_2d_get_x1_delta(this%spline) !this%spline%x1_delta
361  delta_eta2 = sll_f_cubic_spline_2d_get_x2_delta(this%spline) !this%spline%x2_delta
362 
363  call sll_s_cubic_spline_2d_compute_interpolant(data_in, this%spline)
364 
365  if (this%bc_type1 == sll_p_periodic .and. &
366  this%bc_type2 == sll_p_periodic) then
367 
368  do j = 1, num_points2
369  do i = 1, num_points1
370  eta1 = eta1_min + (i - 1)*delta_eta1
371  eta2 = eta2_min + (j - 1)*delta_eta2
372  eta1 = eta1_min + &
373  modulo(eta1 - eta1_min + alpha1(i, j), eta1_max - eta1_min)
374  eta2 = eta2_min + &
375  modulo(eta2 - eta2_min + alpha2(i, j), eta2_max - eta2_min)
376  data_out(i, j) = this%interpolate_from_interpolant_value(eta1, eta2)
377  end do
378  end do
379 
380  else if (this%bc_type1 == sll_p_hermite .and. &
381  this%bc_type2 == sll_p_hermite) then
382 
383  do j = 1, num_points2
384  do i = 1, num_points1
385  eta1 = eta1_min + (i - 1)*delta_eta1 + alpha1(i, j)
386  eta2 = eta2_min + (j - 1)*delta_eta2 + alpha2(i, j)
387  eta1 = min(eta1, eta1_max)
388  eta2 = min(eta2, eta2_max)
389  eta1 = max(eta1, eta1_min)
390  eta2 = max(eta2, eta2_min)
391  data_out(i, j) = this%interpolate_from_interpolant_value(eta1, eta2)
392  end do
393  end do
394 
395  else
396 
397  do j = 1, num_points2
398  do i = 1, num_points1
399  eta1 = eta1_min + (i - 1)*delta_eta1 + alpha1(i, j)
400  eta2 = eta2_min + (j - 1)*delta_eta2 + alpha2(i, j)
401  sll_assert(eta1_min <= eta1 .and. eta1 <= eta1_max)
402  sll_assert(eta2_min <= eta2 .and. eta2 <= eta2_max)
403  data_out(i, j) = this%interpolate_from_interpolant_value(eta1, eta2)
404  end do
405  end do
406  end if
407  end subroutine spline_interpolate2d_disp
408 
409  subroutine set_coefficients_cs2d( &
410  interpolator, &
411  coeffs_1d, &
412  coeffs_2d, &
413  coeff2d_size1, &
414  coeff2d_size2, &
415  knots1, &
416  size_knots1, &
417  knots2, &
418  size_knots2)
419  class(sll_t_cubic_spline_interpolator_2d), intent(inout) :: interpolator
420  sll_real64, dimension(:), intent(in), optional :: coeffs_1d
421  sll_real64, dimension(:, :), intent(in), optional :: coeffs_2d
422  ! size coeffs 2D
423  sll_int32, intent(in), optional :: coeff2d_size1
424  sll_int32, intent(in), optional :: coeff2d_size2
425  sll_real64, dimension(:), intent(in), optional :: knots1
426  sll_real64, dimension(:), intent(in), optional :: knots2
427  sll_int32, intent(in), optional :: size_knots1
428  sll_int32, intent(in), optional :: size_knots2
429  print *, 'set_coefficients_cs2d(): ERROR: This function has not been ', &
430  'implemented yet.'
431  print *, interpolator%npts1
432  if (present(coeffs_1d)) then
433  print *, 'coeffs_1d present but not used'
434  end if
435  if (present(coeffs_2d)) then
436  print *, 'coeffs_2d present but not used'
437  end if
438  if (present(coeff2d_size1)) then
439  print *, 'coeff2d_size1 present but not used'
440  end if
441  if (present(coeff2d_size2)) then
442  print *, 'coeff2d_size2 present but not used'
443  end if
444  if (present(knots1)) then
445  print *, 'knots1 present but not used'
446  end if
447  if (present(knots2)) then
448  print *, 'knots2 present but not used'
449  end if
450  if (present(size_knots1)) then
451  print *, 'size_knots1 present but not used'
452  end if
453  if (present(size_knots2)) then
454  print *, 'size_knots2 present but not used'
455  end if
456 
457  stop
458  end subroutine !set_coefficients_cs2d
459 
460 !!$ subroutine compute_spl_coeff_cs2d(interpolator, &
461 !!$ data_array, &
462 !!$ eta1_coords, &
463 !!$ size_eta1_coords, &
464 !!$ eta2_coords, &
465 !!$ size_eta2_coords )
466 !!$ class(sll_t_cubic_spline_interpolator_2d), intent(inout) :: interpolator
467 !!$ sll_real64, dimension(:,:), intent(in) :: data_array
468 !!$ sll_real64, dimension(:), intent(in),optional :: eta1_coords
469 !!$ sll_real64, dimension(:), intent(in),optional :: eta2_coords
470 !!$ sll_int32, intent(in), optional :: size_eta1_coords
471 !!$ sll_int32, intent(in),optional :: size_eta2_coords
472 !!$
473 !!$ print *, 'compute_coefficients_cs2d(): ERROR: This function has not been',&
474 !!$ 'implemented yet.'
475 !!$ stop
476 !!$ end subroutine compute_spl_coeff_cs2d
477 
478  function get_coefficients_cs2d(interpolator)
479  class(sll_t_cubic_spline_interpolator_2d), intent(in) :: interpolator
480  sll_real64, dimension(:, :), pointer :: get_coefficients_cs2d
481 
482  print *, 'get_coefficients_cs2d(): ERROR: This function has not been ', &
483  'implemented yet.'
484  get_coefficients_cs2d => null()
485  print *, interpolator%npts1
486  stop
487  end function get_coefficients_cs2d
488 
489  function coefficients_are_set_cs2d(interpolator) result(res)
490  class(sll_t_cubic_spline_interpolator_2d), intent(in) :: interpolator
491  logical :: res
492  res = .false.
493  print *, 'coefficients_are_set_cs2d(): this function has not been implemented yet.'
494  print *, '#', interpolator%npts1
495  !stop
496  end function coefficients_are_set_cs2d
497 
Class for the cubic spline sll_c_interpolator_2d.
subroutine compute_interpolants_cs2d(interpolator, data_array, eta1_coords, size_eta1_coords, eta2_coords, size_eta2_coords)
subroutine spline_interpolate2d(this, num_points1, num_points2, data_in, eta1, eta2, data_out)
subroutine initialize_cs2d_interpolator(interpolator, npts1, npts2, eta1_min, eta1_max, eta2_min, eta2_max, eta1_bc_type, eta2_bc_type, const_eta1_min_slope, const_eta1_max_slope, const_eta2_min_slope, const_eta2_max_slope, eta1_min_slopes, eta1_max_slopes, eta2_min_slopes, eta2_max_slopes)
function interpolate_deriv2_cs2d(interpolator, eta1, eta2)
type(sll_t_cubic_spline_interpolator_2d) function, pointer, public sll_f_new_cubic_spline_interpolator_2d(npts1, npts2, eta1_min, eta1_max, eta2_min, eta2_max, eta1_bc_type, eta2_bc_type, const_eta1_min_slope, const_eta1_max_slope, const_eta2_min_slope, const_eta2_max_slope, eta1_min_slopes, eta1_max_slopes, eta2_min_slopes, eta2_max_slopes)
Function that return a pointer to a cubic spline interpolator 2d object. The result can be the target...
subroutine set_coefficients_cs2d(interpolator, coeffs_1d, coeffs_2d, coeff2d_size1, coeff2d_size2, knots1, size_knots1, knots2, size_knots2)
function interpolate_deriv1_cs2d(interpolator, eta1, eta2)
subroutine spline_interpolate2d_disp(this, num_points1, num_points2, data_in, alpha1, alpha2, data_out)
function interpolate_value_cs2d(interpolator, eta1, eta2)
Provides capabilities for data and derivative interpolation with cubic B-splines and different bounda...
subroutine, public sll_s_cubic_spline_2d_compute_interpolant(data, spline)
Computes the spline coefficients for the given data. The coeffcients are first computed in the second...
real(kind=f64) function, public sll_f_cubic_spline_2d_eval_deriv_x2(x1, x2, spline)
Returns the interpolated value of the derivative.
real(kind=f64) function, public sll_f_cubic_spline_2d_eval_deriv_x1(x1, x2, spline)
Returns the interpolated value of the derivative in the x1 direction at the point (x1,...
real(kind=f64) function, public sll_f_cubic_spline_2d_eval(x1, x2, spline)
Returns the interpolated value of the image of the point (x1,x2) using the spline decomposition store...
subroutine, public sll_s_cubic_spline_2d_init(this, num_pts_x1, num_pts_x2, x1_min, x1_max, x2_min, x2_max, x1_bc_type, x2_bc_type, const_slope_x1_min, const_slope_x1_max, const_slope_x2_min, const_slope_x2_max, x1_min_slopes, x1_max_slopes, x2_min_slopes, x2_max_slopes)
subroutine, public sll_s_cubic_spline_2d_free(spline)
abstract data type for 2d interpolation
The spline-based interpolator is only a wrapper around the capabilities of the cubic splines.
Basic type for one-dimensional cubic spline data.
Base class/basic interface for 2D interpolators.
    Report Typos and Errors