Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_fcisl_toroidal.F90
Go to the documentation of this file.
1 module sll_m_fcisl_toroidal
2 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 #include "sll_errors.h"
4 #include "sll_memory.h"
5 #include "sll_working_precision.h"
6 
7  use sll_m_common_array_initializers, only: &
8  sll_i_scalar_initializer_2d
9 
10  use sll_m_constants, only: &
11  sll_p_pi
12 
13  use sll_m_hermite_interpolation_2d, only: &
14  sll_s_compute_hermite_derivatives_periodic1, &
15  sll_s_compute_w_hermite, &
16  sll_s_localize_per
17 
18  use sll_m_lagrange_interpolation, only: &
19  sll_f_lagrange_interpolate
20 
21  implicit none
22 
23  public :: &
24  sll_s_compute_analytic_field, &
25  sll_s_compute_euler_field, &
26  sll_s_compute_feet_analytic, &
27  sll_s_compute_feet_euler, &
28  sll_s_compute_feet_rk4, &
29  sll_f_compute_inverse_invr_integral, &
30  sll_f_compute_invr_integral, &
31  sll_s_compute_linspace, &
32  sll_s_compute_modulo_vect, &
33  sll_s_compute_modulo_vect2d_inplace, &
34  sll_s_compute_rk4_field, &
35  sll_s_compute_time_points, &
36  sll_f_interpolate2d_toroidal
37 
38  private
39 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
40 
41 !we consider i this test a case where magnetic filed
42 !is not straght lines in theta phi
43 !there is a change
44 !however everything can be analitically computed
45 !we write here some routines
46 !that are useful for making the test
47 !which is test_fcisl_toroidal
48 
49 contains
50 
51  subroutine sll_s_compute_time_points(dt, num_time_points, time_points)
52  sll_real64, intent(in) :: dt
53  sll_int32, intent(in) :: num_time_points
54  sll_real64, dimension(:), intent(out) :: time_points
55 
56  sll_int32 :: i
57 
58  if (size(time_points) < num_time_points) then
59  print *, '#bad size for time_points'
60  print *, '#at line/file:', __line__, __file__
61  end if
62 
63  do i = 1, num_time_points
64  time_points(i) = real(i - 1, f64)*dt
65  end do
66 
67  end subroutine sll_s_compute_time_points
68 
69  subroutine sll_s_compute_euler_field( &
70  R0, &
71  time_points, &
72  num_time_points, &
73  psipr, &
74  F0, &
75  smallr, &
76  theta, &
77  phi, &
78  theta0, &
79  phi0)
80  sll_real64, intent(in) :: r0
81  sll_int32, intent(in) :: num_time_points
82  sll_real64, intent(in) :: psipr
83  sll_real64, intent(in) :: f0
84  sll_real64, intent(in) :: smallr
85  sll_real64, dimension(:), intent(in) :: time_points
86  sll_real64, dimension(:), intent(out) :: theta
87  sll_real64, dimension(:), intent(out) :: phi
88  sll_real64, intent(in) :: theta0
89  sll_real64, intent(in) :: phi0
90  sll_real64 :: dt
91 
92  sll_real64 :: bigr
93  !sll_real64 :: ph
94  !sll_real64 :: th
95  sll_int32 :: i
96 
97  theta(1) = theta0
98  phi(1) = phi0
99  do i = 1, num_time_points - 1
100  dt = time_points(i + 1) - time_points(i)
101  bigr = r0 + smallr*cos(theta(i))
102  theta(i + 1) = theta(i) - psipr/smallr*dt
103  phi(i + 1) = phi(i) + f0/bigr*dt
104  end do
105 
106  end subroutine sll_s_compute_euler_field
107 
108  subroutine sll_s_compute_rk4_field( &
109  R0, &
110  time_points, &
111  num_time_points, &
112  psipr, &
113  F0, &
114  smallr, &
115  theta, &
116  phi, &
117  theta0, &
118  phi0)
119 
120  sll_real64, intent(in) :: r0
121  sll_int32, intent(in) :: num_time_points
122  sll_real64, intent(in) :: psipr
123  sll_real64, intent(in) :: f0
124  sll_real64, intent(in) :: smallr
125  sll_real64, dimension(:), intent(in) :: time_points
126  sll_real64, dimension(:), intent(out) :: theta
127  sll_real64, dimension(:), intent(out) :: phi
128  sll_real64, intent(in) :: theta0
129  sll_real64, intent(in) :: phi0
130  procedure(sll_i_scalar_initializer_2d), pointer :: func1
131  sll_real64 :: params1(2)
132  procedure(sll_i_scalar_initializer_2d), pointer :: func2
133  sll_real64 :: params2(3)
134  sll_real64 :: dt
135  !sll_real64 :: bigr
136  !sll_real64 :: ph
137  !sll_real64 :: th
138  sll_int32 :: i
139  sll_real64 :: k1(4)
140  sll_real64 :: k2(4)
141 
142  func1 => toroidal_1
143  params1(1) = smallr
144  params1(2) = psipr
145  func2 => toroidal_2
146  params2(1) = smallr
147  params2(2) = r0
148  params2(3) = f0
149  theta(1) = theta0
150  phi(1) = phi0
151  do i = 1, num_time_points - 1
152  dt = time_points(i + 1) - time_points(i)
153  k1(1) = func1(theta(i), phi(i), params1)
154  k2(1) = func2(theta(i), phi(i), params2)
155  k1(2) = func1(theta(i) + 0.5_f64*dt*k1(1), phi(i) + 0.5_f64*dt*k2(1), params1)
156  k2(2) = func2(theta(i) + 0.5_f64*dt*k1(1), phi(i) + 0.5_f64*dt*k2(1), params2)
157  k1(3) = func1(theta(i) + 0.5_f64*dt*k1(2), phi(i) + 0.5_f64*dt*k2(2), params1)
158  k2(3) = func2(theta(i) + 0.5_f64*dt*k1(2), phi(i) + 0.5_f64*dt*k2(2), params2)
159  k1(4) = func1(theta(i) + dt*k1(3), phi(i) + dt*k2(3), params1)
160  k2(4) = func2(theta(i) + dt*k1(3), phi(i) + dt*k2(3), params2)
161  theta(i + 1) = theta(i) + (dt/6._f64)*(k1(1) + 2._f64*(k1(2) + k1(3)) + k1(4))
162  phi(i + 1) = phi(i) + (dt/6._f64)*(k2(1) + 2._f64*(k2(2) + k2(3)) + k2(4))
163  end do
164 
165  end subroutine sll_s_compute_rk4_field
166 
167  function toroidal_1(theta, phi, params) result(res)
168  sll_real64 :: res
169  sll_real64, intent(in) :: theta
170  sll_real64, intent(in) :: phi
171  sll_real64, dimension(:), intent(in) :: params
172 
173  sll_real64 :: psipr
174  sll_real64 :: smallr
175 
176  smallr = params(1)! default: 1.0_f64
177  psipr = params(2) ! default: 1.0_f64
178 
179  res = -psipr/smallr
180  return !PN ADD TO PREVENT WARNING
181  print *, phi, theta
182  end function toroidal_1
183 
184  function toroidal_2(theta, phi, params) result(res)
185  sll_real64 :: res
186  sll_real64, intent(in) :: theta
187  sll_real64, intent(in) :: phi
188  sll_real64, dimension(:), intent(in) :: params
189 
190  sll_real64 :: bigr
191  sll_real64 :: f0
192  sll_real64 :: r0
193  sll_real64 :: smallr
194 
195  smallr = params(1) ! default: 1.0_f64
196  r0 = params(2) ! default: 1.0_f64
197  f0 = params(3) ! default: 1.0_f64
198 
199  bigr = r0 + smallr*cos(theta)
200  res = f0/bigr
201 
202  return !PN ADD TO PREVENT WARNING
203  print *, phi
204  end function toroidal_2
205 
206  subroutine sll_s_compute_analytic_field( &
207  R0, &
208  time_points, &
209  num_time_points, &
210  psipr, &
211  F0, &
212  smallr, &
213  theta, &
214  phi, &
215  theta0, &
216  phi0)
217  sll_real64, intent(in) :: r0
218  sll_int32, intent(in) :: num_time_points
219  sll_real64, intent(in) :: psipr
220  sll_real64, intent(in) :: f0
221  sll_real64, intent(in) :: smallr
222  sll_real64, dimension(:), intent(in) :: time_points
223  sll_real64, dimension(:), intent(out) :: theta
224  sll_real64, dimension(:), intent(out) :: phi
225  sll_real64, intent(in) :: theta0
226  sll_real64, intent(in) :: phi0
227  !sll_real64 :: dt
228 
229  !sll_real64 :: bigr
230  !sll_real64 :: ph
231  !sll_real64 :: th
232  sll_int32 :: i
233 
234  do i = 1, num_time_points
235  theta(i) = theta0 - psipr/smallr*(time_points(i) - time_points(1))
236  phi(i) = phi0 -f0*smallr/psipr*(sll_f_compute_invr_integral(r0,smallr,theta(i))-sll_f_compute_invr_integral(r0,smallr,theta0))
237  end do
238 
239  end subroutine sll_s_compute_analytic_field
240 
241  function sll_f_compute_invr_integral(R0, smallr, x) result(res)
242  sll_real64 :: res
243  sll_real64, intent(in) :: r0
244  sll_real64, intent(in) :: smallr
245  sll_real64, intent(in) :: x
246 
247  sll_real64 :: xx
248  sll_real64 :: k
249  sll_real64 :: alpha
250  !computes int(1/(R0+smallr*cos(u)),u=0..x)
251  !we suppose that R0>r
252 
253  xx = 0.5_f64 + x/(2._f64*sll_p_pi)
254  k = real(floor(xx), f64)
255  alpha = xx - k
256  if (alpha < 1.e-12) then
257  res = -0.5_f64*sll_p_pi
258  else
259  res = tan((-0.5_f64 + alpha)*sll_p_pi)
260  res = atan((r0 - smallr)/sqrt(r0**2 - smallr**2)*res)
261  end if
262  res = res + k*sll_p_pi
263  res = 2._f64/(sqrt(r0**2 - smallr**2))*res
264 
265  end function sll_f_compute_invr_integral
266 
267  function sll_f_compute_inverse_invr_integral(R0, smallr, x) result(res)
268  sll_real64 :: res
269  sll_real64, intent(in) :: r0
270  sll_real64, intent(in) :: smallr
271  sll_real64, intent(in) :: x
272 
273  sll_real64 :: xx
274  sll_real64 :: k
275  sll_real64 :: alpha
276  !computes inverse of function int(1/(R0+smallr*cos(u)),u=0..x)
277  !we suppose that R0>r
278 
279  xx = 0.5_f64*x*sqrt(r0**2 - smallr**2)
280  xx = 0.5_f64 + xx/sll_p_pi
281  k = real(floor(xx), f64)
282  alpha = xx - k
283  if (alpha < 1.e-12) then
284  res = -0.5_f64*sll_p_pi
285  else
286  res = tan((-0.5_f64 + alpha)*sll_p_pi)
287  res = atan(sqrt(r0**2 - smallr**2)/(r0 - smallr)*res)
288  end if
289  res = res + k*sll_p_pi
290  res = 2._f64*res
291 
292  end function sll_f_compute_inverse_invr_integral
293 
294  subroutine compute_modulo(x, L)
295  sll_real64, intent(inout) :: x
296  sll_real64, intent(in) :: l
297 
298  x = x/l + 1._f64
299  x = x - real(floor(x), f64)
300  x = l*x
301 
302  end subroutine compute_modulo
303 
304  subroutine sll_s_compute_modulo_vect(in, out, num_points, L)
305  sll_real64, dimension(:), intent(in) :: in
306  sll_real64, dimension(:), intent(out) :: out
307  sll_int32, intent(in) :: num_points
308  sll_real64, intent(in) :: l
309 
310  sll_int32 :: i
311 
312  if (size(in) < num_points) then
313  print *, '#problem of size of in', size(in), num_points
314  print *, '#at line/file', __line__, __file__
315  stop
316  end if
317  if (size(out) < num_points) then
318  print *, '#problem of size of out', size(out), num_points
319  print *, '#at line/file', __line__, __file__
320  stop
321  end if
322 
323  do i = 1, num_points
324  out(i) = in(i)
325  call compute_modulo(out(i), l)
326  end do
327 
328  end subroutine sll_s_compute_modulo_vect
329 
330  subroutine sll_s_compute_modulo_vect2d_inplace(inout, num_points1, num_points2, L)
331  sll_real64, dimension(:, :), intent(inout) :: inout
332  sll_int32, intent(in) :: num_points1
333  sll_int32, intent(in) :: num_points2
334  sll_real64, intent(in) :: l
335 
336  sll_int32 :: i1
337  sll_int32 :: i2
338 
339  if (size(inout, 1) < num_points1) then
340  print *, '#problem of size1 of inout', size(inout, 1), num_points1
341  print *, '#at line/file', __line__, __file__
342  stop
343  end if
344  if (size(inout, 2) < num_points2) then
345  print *, '#problem of size2 of inout', size(inout, 2), num_points2
346  print *, '#at line/file', __line__, __file__
347  stop
348  end if
349 
350  do i2 = 1, num_points2
351  do i1 = 1, num_points1
352  call compute_modulo(inout(i1, i2), l)
353  end do
354  end do
355 
356  end subroutine sll_s_compute_modulo_vect2d_inplace
357 
358  subroutine sll_s_compute_linspace(array, xmin, xmax, Npts)
359  sll_real64, dimension(:), intent(out) :: array
360  sll_real64, intent(in) :: xmin
361  sll_real64, intent(in) :: xmax
362  sll_int32, intent(in) :: npts
363 
364  sll_int32 :: i
365  sll_real64 :: dx
366 
367  dx = (xmax - xmin)/real(npts - 1, f64)
368 
369  do i = 1, npts
370  array(i) = xmin + real(i - 1, f64)*dx
371  end do
372 
373  end subroutine sll_s_compute_linspace
374 
375  subroutine sll_s_compute_feet_euler(theta_in, phi_in, num_points1, num_points2, dt, params, theta_out, phi_out)
376  sll_real64, dimension(:), intent(in) :: theta_in
377  sll_real64, dimension(:), intent(in) :: phi_in
378  sll_int32, intent(in) :: num_points1
379  sll_int32, intent(in) :: num_points2
380  sll_real64, intent(in) :: dt
381  sll_real64, dimension(:), intent(in) :: params
382  sll_real64, dimension(:, :), intent(out) :: theta_out
383  sll_real64, dimension(:, :), intent(out) :: phi_out
384 
385  sll_real64 :: r0
386  sll_real64 :: psipr
387  sll_real64 :: f0
388  sll_real64 :: smallr
389 
390  sll_real64 :: bigr
391  !sll_real64 :: ph
392  !sll_real64 :: th
393  sll_int32 :: i
394  sll_int32 :: j
395 
396  if (size(params) < 4) then
397  print *, '#size of params not good', size(params)
398  print *, '#at line/file', __line__, __file__
399  stop
400  end if
401 
402  r0 = params(1)
403  psipr = params(2)
404  f0 = params(3)
405  smallr = params(4)
406 
407  do j = 1, num_points2
408  do i = 1, num_points1
409  bigr = r0 + smallr*cos(theta_in(i))
410  theta_out(i, j) = theta_in(i) - psipr/smallr*dt
411  phi_out(i, j) = phi_in(j) + f0/bigr*dt
412  end do
413  end do
414 
415  end subroutine sll_s_compute_feet_euler
416 
417  subroutine sll_s_compute_feet_rk4(theta_in, phi_in, num_points1, num_points2, dt, params, theta_out, phi_out)
418  sll_real64, dimension(:), intent(in) :: theta_in
419  sll_real64, dimension(:), intent(in) :: phi_in
420  sll_int32, intent(in) :: num_points1
421  sll_int32, intent(in) :: num_points2
422  sll_real64, intent(in) :: dt
423  sll_real64, dimension(:), intent(in) :: params
424  sll_real64, dimension(:, :), intent(out) :: theta_out
425  sll_real64, dimension(:, :), intent(out) :: phi_out
426 
427  sll_real64 :: r0
428  sll_real64 :: psipr
429  sll_real64 :: f0
430  sll_real64 :: smallr
431 
432  !sll_real64 :: bigr
433  !sll_real64 :: ph
434  !sll_real64 :: th
435  sll_int32 :: i
436  sll_int32 :: j
437 
438  procedure(sll_i_scalar_initializer_2d), pointer :: func1
439  sll_real64 :: params1(2)
440  procedure(sll_i_scalar_initializer_2d), pointer :: func2
441  sll_real64 :: params2(3)
442  sll_real64 :: k1(4)
443  sll_real64 :: k2(4)
444 
445  if (size(params) < 4) then
446  print *, '#size of params not good', size(params)
447  print *, '#at line/file', __line__, __file__
448  stop
449  end if
450 
451  r0 = params(1)
452  psipr = params(2)
453  f0 = params(3)
454  smallr = params(4)
455 
456  func1 => toroidal_1
457  params1(1) = smallr
458  params1(2) = psipr
459  func2 => toroidal_2
460  params2(1) = smallr
461  params2(2) = r0
462  params2(3) = f0
463 
464  do j = 1, num_points2
465  do i = 1, num_points1
466  k1(1) = func1(theta_in(i), phi_in(j), params1)
467  k2(1) = func2(theta_in(i), phi_in(j), params2)
468  k1(2) = func1(theta_in(i) + 0.5_f64*dt*k1(1), phi_in(j) + 0.5_f64*dt*k2(1), params1)
469  k2(2) = func2(theta_in(i) + 0.5_f64*dt*k1(1), phi_in(j) + 0.5_f64*dt*k2(1), params2)
470  k1(3) = func1(theta_in(i) + 0.5_f64*dt*k1(2), phi_in(j) + 0.5_f64*dt*k2(2), params1)
471  k2(3) = func2(theta_in(i) + 0.5_f64*dt*k1(2), phi_in(j) + 0.5_f64*dt*k2(2), params2)
472  k1(4) = func1(theta_in(i) + dt*k1(3), phi_in(j) + dt*k2(3), params1)
473  k2(4) = func2(theta_in(i) + dt*k1(3), phi_in(j) + dt*k2(3), params2)
474  theta_out(i, j) = theta_in(i) + (dt/6._f64)*(k1(1) + 2._f64*(k1(2) + k1(3)) + k1(4))
475  phi_out(i, j) = phi_in(j) + (dt/6._f64)*(k2(1) + 2._f64*(k2(2) + k2(3)) + k2(4))
476  end do
477  end do
478 
479  end subroutine sll_s_compute_feet_rk4
480 
481  subroutine sll_s_compute_feet_analytic(theta_in, phi_in, num_points1, num_points2, dt, params, theta_out, phi_out)
482  sll_real64, dimension(:), intent(in) :: theta_in
483  sll_real64, dimension(:), intent(in) :: phi_in
484  sll_int32, intent(in) :: num_points1
485  sll_int32, intent(in) :: num_points2
486  sll_real64, intent(in) :: dt
487  sll_real64, dimension(:), intent(in) :: params
488  sll_real64, dimension(:, :), intent(out) :: theta_out
489  sll_real64, dimension(:, :), intent(out) :: phi_out
490 
491  sll_real64 :: r0
492  sll_real64 :: psipr
493  sll_real64 :: f0
494  sll_real64 :: smallr
495 
496  sll_real64 :: bigr
497  !sll_real64 :: ph
498  !sll_real64 :: th
499  sll_int32 :: i
500  sll_int32 :: j
501 
502  if (size(params) < 4) then
503  print *, '#size of params not good', size(params)
504  print *, '#at line/file', __line__, __file__
505  stop
506  end if
507 
508  r0 = params(1)
509  psipr = params(2)
510  f0 = params(3)
511  smallr = params(4)
512 
513  do j = 1, num_points2
514  do i = 1, num_points1
515  bigr = r0 + smallr*cos(theta_in(i))
516  theta_out(i, j) = theta_in(i) - psipr/smallr*dt
517  phi_out(i, j) = phi_in(j) &
518  - f0*smallr/psipr*( &
519  sll_f_compute_invr_integral(r0, smallr, theta_out(i, j)) &
520  - sll_f_compute_invr_integral(r0, smallr, theta_in(i)))
521  end do
522  end do
523 
524  end subroutine sll_s_compute_feet_analytic
525 
526  function sll_f_interpolate2d_toroidal( &
527  num_points1, &
528  num_points2, &
529  data_in, &
530  eta1, &
531  eta2, &
532  params) &
533  result(data_out)
534  sll_int32, intent(in) :: num_points1
535  sll_int32, intent(in) :: num_points2
536  sll_real64, dimension(:, :), intent(in) :: eta1
537  sll_real64, dimension(:, :), intent(in) :: eta2
538  sll_real64, dimension(:, :), intent(in) :: data_in
539  sll_real64, dimension(num_points1, num_points2) :: data_out
540  sll_real64, dimension(:), intent(in) :: params
541 
542  sll_real64 :: r0
543  sll_real64 :: psipr
544  sll_real64 :: f0
545  sll_real64 :: smallr
546  sll_int32 :: hermite_p
547  sll_int32 :: lag_r
548  sll_int32 :: lag_s
549  sll_int32 :: lag_p
550  sll_int32 :: ierr
551  sll_real64, dimension(:, :, :), allocatable :: hermite_buf
552  sll_int32 :: i
553  sll_int32 :: j
554  sll_real64 :: eta1_min
555  sll_real64 :: eta1_max
556  sll_real64 :: eta2_min
557  sll_real64 :: eta2_max
558  !sll_real64 :: eta1_star
559  !sll_real64 :: eta2_star
560  sll_real64, dimension(:), allocatable :: lag_buf
561  sll_real64, dimension(:), allocatable :: lag_x
562  sll_real64 :: eta1_loc
563  sll_real64 :: eta2_loc
564  sll_int32 :: i0
565  sll_int32 :: j0
566  sll_int32 :: j0_loc
567  sll_int32 :: jj
568  sll_real64 :: w(4)
569  sll_real64 :: tmp1
570  sll_real64 :: tmp2
571  !we define aligned interpolation
572  !Lagrange interpolation in aligned direction (lag_r,lag_s)
573  !Hermite interpolation in theta direction (hermite_p)
574  !alignement with respect to magnetic field lines
575  !given analytically
576 
577  r0 = params(1)
578  psipr = params(2)
579  f0 = params(3)
580  smallr = params(4)
581  hermite_p = int(params(5), i32)
582  lag_r = int(params(6), i32)
583  lag_s = int(params(7), i32)
584  eta1_min = params(8)
585  eta1_max = params(9)
586  eta2_min = params(10)
587  eta2_max = params(11)
588 
589  lag_p = lag_s - lag_r
590 
591  sll_allocate(hermite_buf(3, num_points1, num_points2), ierr)
592  call sll_s_compute_hermite_derivatives_periodic1( &
593  data_in, &
594  num_points1, &
595  num_points2, &
596  hermite_p, &
597  hermite_buf)
598 
599  !print *,'hermite_buf=',maxval(abs(hermite_buf(1,:,:)-1._f64)), &
600  ! maxval(abs(hermite_buf(2,:,:))),maxval(abs(hermite_buf(3,:,:)))
601  !stop
602 
603  sll_allocate(lag_buf(lag_r:lag_s), ierr)
604  sll_allocate(lag_x(lag_r:lag_s), ierr)
605  do i = lag_r, lag_s
606  lag_x(i) = real(i, f64)
607  end do
608 
609  tmp1 = (eta2_max - eta2_min)/real(num_points2 - 1, f64)*(-psipr)/(f0*smallr)
610  do j = 1, num_points2
611  do i = 1, num_points1
612  eta2_loc = eta2(i, j)
613  call sll_s_localize_per(j0, eta2_loc, eta2_min, eta2_max, num_points2 - 1)
614  tmp2 = sll_f_compute_invr_integral(r0, smallr, eta1(i, j))
615  tmp2 = tmp2 - eta2_loc*tmp1
616  do jj = lag_r, lag_s
617  j0_loc = modulo(j0 + jj + num_points2 - 1, num_points2 - 1) + 1
618  eta1_loc = sll_f_compute_inverse_invr_integral(r0, smallr, real(jj, f64)*tmp1 + tmp2)
619  !phi(t) = phi(0)-F0r/psipr*(L(th(t))-L(th(0)))
620  !Hermite interpolation
621  call sll_s_localize_per(i0, eta1_loc, eta1_min, eta1_max, num_points1 - 1)
622  i0 = i0 + 1
623  w(1) = (2._f64*eta1_loc + 1)*(1._f64 - eta1_loc)*(1._f64 - eta1_loc)
624  w(2) = eta1_loc*eta1_loc*(3._f64 - 2._f64*eta1_loc)
625  w(3) = eta1_loc*(1._f64 - eta1_loc)*(1._f64 - eta1_loc)
626  w(4) = eta1_loc*eta1_loc*(eta1_loc - 1._f64)
627 
628  !if(abs(eta2_loc)<1.e-12)then
629  ! print *,eta1_loc,eta2_loc,j0,jj
630  !endif
631  lag_buf(jj) = w(1)*hermite_buf(1, i0, j0_loc) &
632  + w(2)*hermite_buf(1, i0 + 1, j0_loc) &
633  + w(3)*hermite_buf(2, i0, j0_loc) &
634  + w(4)*hermite_buf(3, i0, j0_loc)
635  !if(abs(lag_buf(jj)-1._f64)>1.e-12)then
636  ! print *,i0,j0_loc,lag_buf(jj)-1._f64,hermite_buf(1,i0,j0_loc),hermite_buf(2,i0,j0_loc),hermite_buf(3,i0,j0_loc),&
637  ! hermite_buf(1,i0+1,j0_loc)
638  ! stop
639  !endif
640  end do
641 
642  data_out(i, j) = sll_f_lagrange_interpolate( &
643  eta2_loc, &
644  lag_p, &
645  lag_x(lag_r:lag_s), &
646  lag_buf(lag_r:lag_s))
647  end do
648  end do
649  end function sll_f_interpolate2d_toroidal
650 
651  subroutine compute_w_hermite_aligned( &
652  w, &
653  w_cell, &
654  num_points1, &
655  r, &
656  s, &
657  eta1_pos, &
658  eta1_min, &
659  eta1_max)
660  sll_real64, dimension(:, :, :), intent(out) :: w
661  sll_int32, dimension(:, :), intent(out) :: w_cell
662  sll_int32, intent(in) :: num_points1
663  sll_int32, intent(in) :: r
664  sll_int32, intent(in) :: s
665  sll_real64, dimension(:, :), intent(in) :: eta1_pos
666  sll_real64, intent(in) :: eta1_min
667  sll_real64, intent(in) :: eta1_max
668  sll_int32 :: i
669  !sll_int32 :: j
670  sll_int32 :: ell
671  sll_real64 :: w_loc_tmp(r:s)
672  sll_real64 :: w_loc(s - r)
673  sll_int32 :: i0
674  sll_real64 :: eta1_loc
675  sll_real64 :: w_basis(4)
676  sll_int32 :: ii
677  character(len=*), parameter :: fun = "compute_w_hermite_aligned"
678 
679  if (r > 0 .or. s < 0) then
680  print *, '#r,s=', r, s
681  print *, '#not treated'
682  sll_error(fun, '#bad value of r and s')
683  end if
684  if (size(w, 1) < 4) then
685  sll_error(fun, '#bad size1 for w')
686  end if
687  if (size(w, 2) < s - r) then
688  sll_error(fun, '#bad size2 for w')
689  end if
690  if (size(w, 3) < num_points1) then
691  sll_error(fun, '#bad size3 for w')
692  end if
693  if (size(eta1_pos, 1) < s - r) then
694  sll_error(fun, '#bad size1 for eta1_pos')
695  end if
696  if (size(eta1_pos, 2) < num_points1) then
697  sll_error(fun, '#bad size2 for eta1_pos')
698  end if
699 
700  call sll_s_compute_w_hermite(w_loc_tmp, r, s)
701  w_loc(1:-r) = w_loc_tmp(r:-1)
702  w_loc(-r + 1:s) = w_loc_tmp(1:s)
703 
704  do i = 1, num_points1
705  do ell = 1, r - s
706  eta1_loc = eta1_pos(ell, i)
707  call sll_s_localize_per(i0, eta1_loc, eta1_min, eta1_max, num_points1 - 1)
708  i0 = i0 + 1
709  w_cell(ell, i) = i0
710  w_basis(1) = (2._f64*eta1_loc + 1)*(1._f64 - eta1_loc)*(1._f64 - eta1_loc)
711  w_basis(2) = eta1_loc*eta1_loc*(3._f64 - 2._f64*eta1_loc)
712  w_basis(3) = eta1_loc*(1._f64 - eta1_loc)*(1._f64 - eta1_loc)
713  w_basis(4) = eta1_loc*eta1_loc*(eta1_loc - 1._f64)
714  do ii = 1, 4
715  w(ii, ell, i) = w_loc(ell)*w_basis(ii)
716  end do
717  end do
718  end do
719 
720  end subroutine compute_w_hermite_aligned
721 
722  subroutine compute_hermite_derivatives_aligned( &
723  f, &
724  num_points1, &
725  num_points2, &
726  p1, &
727  w_aligned_left, &
728  w_cell_aligned_left, &
729  w_aligned_right, &
730  w_cell_aligned_right, &
731  buf)
732  sll_real64, dimension(:, :), intent(in) :: f
733  sll_int32, intent(in) :: num_points1
734  sll_int32, intent(in) :: num_points2
735  sll_int32, intent(in) :: p1
736  sll_real64, dimension(:, :, :), intent(in) :: w_aligned_left
737  sll_real64, dimension(:, :, :), intent(in) :: w_aligned_right
738  sll_int32, dimension(:, :), intent(in) :: w_cell_aligned_left
739  sll_int32, dimension(:, :), intent(in) :: w_cell_aligned_right
740  sll_real64, dimension(:, :, :), intent(out) :: buf
741 
742  sll_real64 :: w_left(-p1/2:(p1 + 1)/2)
743  sll_real64 :: w_right((-p1 + 1)/2:p1/2 + 1)
744  sll_int32 :: r_left
745  sll_int32 :: s_left
746  sll_int32 :: r_right
747  sll_int32 :: s_right
748  sll_int32 :: i
749  sll_int32 :: j
750  sll_real64 :: tmp
751  sll_int32 :: ii
752  sll_int32 :: ind
753 
754  r_left = -p1/2
755  s_left = (p1 + 1)/2
756  r_right = (-p1 + 1)/2
757  s_right = p1/2 + 1
758  call sll_s_compute_w_hermite(w_left, r_left, s_left)
759  if (((2*p1/2) - p1) == 0) then
760  w_right(r_right:s_right) = w_left(r_left:s_left)
761  else
762  w_right(r_right:s_right) = -w_left(s_left:r_left:-1)
763  end if
764 
765  !first compute \partial_1f(eta1_i,eta2_j)
766  do j = 1, num_points2
767  do i = 1, num_points1
768  buf(1, i, j) = f(i, j) !f(0)
769  tmp = 0._f64
770  do ii = r_left, s_left
771  ind = modulo(i + ii - 2 + num_points1, num_points1 - 1) + 1
772  tmp = tmp + w_left(ii)*f(ind, j)
773  end do
774  buf(2, i, j) = tmp !df(0)
775  tmp = 0._f64
776  do ii = r_right, s_right
777  ind = modulo(i + ii - 2 + num_points1, num_points1 - 1) + 1
778  tmp = tmp + w_right(ii)*f(ind, j)
779  end do
780  buf(3, i, j) = tmp !df(1)
781  end do
782  end do
783 
784  !then compute \partial_aligned f(eta1_i,eta2_j)
785 
786  return !PN ADD TO PREVENT WARNING
787  print *, w_aligned_left, &
788  w_cell_aligned_left, &
789  w_aligned_right, &
790  w_cell_aligned_right
791 
792  end subroutine compute_hermite_derivatives_aligned
793 
794 end module sll_m_fcisl_toroidal
sll_m_common_array_initializers::sll_i_scalar_initializer_2d
Definition: sll_m_common_array_initializers.F90:89
sll_m_common_array_initializers
Definition: sll_m_common_array_initializers.F90:1
sll_m_constants
Fortran module where set some physical and mathematical constants.
Definition: sll_m_constants.F90:23
sll_m_constants::sll_p_pi
real(kind=f64), parameter, public sll_p_pi
Definition: sll_m_constants.F90:51
sll_m_fcisl_toroidal
Definition: sll_m_fcisl_toroidal.F90:1
sll_m_fcisl_toroidal::sll_s_compute_time_points
subroutine, public sll_s_compute_time_points(dt, num_time_points, time_points)
Definition: sll_m_fcisl_toroidal.F90:52
sll_m_fcisl_toroidal::sll_f_interpolate2d_toroidal
real(kind=f64) function, dimension(num_points1, num_points2), public sll_f_interpolate2d_toroidal(num_points1, num_points2, data_in, eta1, eta2, params)
Definition: sll_m_fcisl_toroidal.F90:534
sll_m_fcisl_toroidal::sll_s_compute_analytic_field
subroutine, public sll_s_compute_analytic_field(R0, time_points, num_time_points, psipr, F0, smallr, theta, phi, theta0, phi0)
Definition: sll_m_fcisl_toroidal.F90:217
sll_m_fcisl_toroidal::sll_s_compute_modulo_vect
subroutine, public sll_s_compute_modulo_vect(in, out, num_points, L)
Definition: sll_m_fcisl_toroidal.F90:305
sll_m_fcisl_toroidal::sll_s_compute_linspace
subroutine, public sll_s_compute_linspace(array, xmin, xmax, Npts)
Definition: sll_m_fcisl_toroidal.F90:359
sll_m_fcisl_toroidal::toroidal_1
real(kind=f64) function toroidal_1(theta, phi, params)
Definition: sll_m_fcisl_toroidal.F90:168
sll_m_fcisl_toroidal::sll_s_compute_feet_rk4
subroutine, public sll_s_compute_feet_rk4(theta_in, phi_in, num_points1, num_points2, dt, params, theta_out, phi_out)
Definition: sll_m_fcisl_toroidal.F90:418
sll_m_fcisl_toroidal::sll_f_compute_inverse_invr_integral
real(kind=f64) function, public sll_f_compute_inverse_invr_integral(R0, smallr, x)
Definition: sll_m_fcisl_toroidal.F90:268
sll_m_fcisl_toroidal::sll_f_compute_invr_integral
real(kind=f64) function, public sll_f_compute_invr_integral(R0, smallr, x)
Definition: sll_m_fcisl_toroidal.F90:242
sll_m_fcisl_toroidal::sll_s_compute_feet_euler
subroutine, public sll_s_compute_feet_euler(theta_in, phi_in, num_points1, num_points2, dt, params, theta_out, phi_out)
Definition: sll_m_fcisl_toroidal.F90:376
sll_m_fcisl_toroidal::sll_s_compute_modulo_vect2d_inplace
subroutine, public sll_s_compute_modulo_vect2d_inplace(inout, num_points1, num_points2, L)
Definition: sll_m_fcisl_toroidal.F90:331
sll_m_fcisl_toroidal::sll_s_compute_rk4_field
subroutine, public sll_s_compute_rk4_field(R0, time_points, num_time_points, psipr, F0, smallr, theta, phi, theta0, phi0)
Definition: sll_m_fcisl_toroidal.F90:119
sll_m_fcisl_toroidal::compute_modulo
subroutine compute_modulo(x, L)
Definition: sll_m_fcisl_toroidal.F90:295
sll_m_fcisl_toroidal::sll_s_compute_euler_field
subroutine, public sll_s_compute_euler_field(R0, time_points, num_time_points, psipr, F0, smallr, theta, phi, theta0, phi0)
Definition: sll_m_fcisl_toroidal.F90:80
sll_m_fcisl_toroidal::compute_hermite_derivatives_aligned
subroutine compute_hermite_derivatives_aligned(f, num_points1, num_points2, p1, w_aligned_left, w_cell_aligned_left, w_aligned_right, w_cell_aligned_right, buf)
Definition: sll_m_fcisl_toroidal.F90:732
sll_m_fcisl_toroidal::sll_s_compute_feet_analytic
subroutine, public sll_s_compute_feet_analytic(theta_in, phi_in, num_points1, num_points2, dt, params, theta_out, phi_out)
Definition: sll_m_fcisl_toroidal.F90:482
sll_m_fcisl_toroidal::compute_w_hermite_aligned
subroutine compute_w_hermite_aligned(w, w_cell, num_points1, r, s, eta1_pos, eta1_min, eta1_max)
Definition: sll_m_fcisl_toroidal.F90:660
sll_m_fcisl_toroidal::toroidal_2
real(kind=f64) function toroidal_2(theta, phi, params)
Definition: sll_m_fcisl_toroidal.F90:185
sll_m_hermite_interpolation_2d
Definition: sll_m_hermite_interpolation_2d.F90:1
sll_m_hermite_interpolation_2d::sll_s_localize_per
subroutine, public sll_s_localize_per(i, x, xmin, xmax, N)
Definition: sll_m_hermite_interpolation_2d.F90:682
sll_m_hermite_interpolation_2d::sll_s_compute_hermite_derivatives_periodic1
subroutine, public sll_s_compute_hermite_derivatives_periodic1(f, num_points1, num_points2, p, buf)
Definition: sll_m_hermite_interpolation_2d.F90:330
sll_m_hermite_interpolation_2d::sll_s_compute_w_hermite
subroutine, public sll_s_compute_w_hermite(w, r, s)
Definition: sll_m_hermite_interpolation_2d.F90:255
sll_m_lagrange_interpolation
Definition: sll_m_lagrange_interpolation.F90:18
sll_m_lagrange_interpolation::sll_f_lagrange_interpolate
real(kind=f64) function, public sll_f_lagrange_interpolate(x, degree, xi, yi)
Definition: sll_m_lagrange_interpolation.F90:56
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:41 for Semi-Lagrangian Library by   doxygen 1.9.1