Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_qn_2d_polar_precompute.F90
Go to the documentation of this file.
1 module sll_m_qn_2d_polar_precompute
2 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 #include "sll_memory.h"
4 #include "sll_working_precision.h"
5 
6 ! use F77_blas, only: &
7 ! dgemm, &
8 ! zgemm
9 
10 ! use F77_fftpack, only: &
11 ! zfftb, &
12 ! zfftf, &
13 ! zffti
14 
15 ! use F77_lapack, only: &
16 ! zgetrf, &
17 ! zgetri
18 
19  use sll_m_constants, only: &
20  sll_p_pi
21 
22  use sll_m_qn_2d_polar, only: &
23  sll_s_compute_splines_coefs_matrix_nat_1d, &
24  sll_s_compute_splines_coefs_matrix_per_1d, &
25  sll_s_contribution_spl, &
26  sll_s_localize_polar, &
27  sll_s_matrix_product_compf, &
28  sll_s_splcoefnat1d0old, &
29  sll_s_splcoefper1d0old
30 
31  implicit none
32 
33  public :: &
34  sll_s_compute_qns_inverse_polar_splines, &
35  sll_s_compute_qns_matrix_polar_splines, &
36  sll_s_solve_qns_polar_splines
37 
38  private
39 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
40 
41 contains
42 
43  subroutine pre_precompute_double_gyroaverage_coeff_polar_splines( &
44  r_min, &
45  r_max, &
46  num_cells_r, &
47  num_cells_theta, &
48  rho, &
49  points, &
50  N_points, &
51  pre_compute_N, &
52  size_pre_compute)
53  sll_real64, intent(in) :: r_min
54  sll_real64, intent(in) :: r_max
55  sll_int32, intent(in) :: num_cells_r
56  sll_int32, intent(in) :: num_cells_theta
57  sll_real64, dimension(:), intent(in) :: rho
58  sll_real64, dimension(:, :), intent(in) :: points
59  sll_int32, intent(in) :: n_points
60  sll_int32, dimension(:), intent(out) :: pre_compute_n
61  sll_int32, intent(out) :: size_pre_compute
62  sll_int32, dimension(:, :), allocatable :: buf
63  sll_int32 ::i, k, ell_1, ell_2, ii(2), s, nb, ind(2)
64  sll_real64::eta_star(2), eta(2), delta_eta(2), x(2)
65  sll_int32 ::error, max_nb
66  sll_real64 :: eta_min(2), eta_max(2)
67  sll_int32 :: nc(2)
68 
69  delta_eta(1) = (r_max - r_min)/real(num_cells_r, f64)
70  delta_eta(2) = 2._f64*sll_p_pi/real(num_cells_theta, f64)
71 
72  nc(1) = num_cells_r
73  nc(2) = num_cells_theta
74 
75  eta_min(1) = r_min
76  eta_min(2) = 0._f64
77  eta_max(1) = r_max
78  eta_max(2) = 2._f64*sll_p_pi
79 
80  sll_allocate(buf(0:num_cells_r + 2, 0:num_cells_theta - 1), error)
81 
82  eta(2) = 0._f64
83 
84  max_nb = 0
85  buf = 0
86  nb = 0
87  do i = 1, nc(1) + 1
88  eta(1) = eta_min(1) + real(i - 1, f64)*delta_eta(1)
89  s = 0
90  do k = 1, n_points
91  x(1) = eta(1)*cos(eta(2)) + rho(i)*points(1, k)
92  x(2) = eta(1)*sin(eta(2)) + rho(i)*points(2, k)
93  call sll_s_localize_polar( &
94  x, &
95  eta_min, &
96  eta_max, &
97  ii, &
98  eta_star, &
99  nc)
100  do ell_2 = -1, 2
101  ind(2) = modulo(ii(2) + ell_2, nc(2))
102  do ell_1 = -1, 2
103  ind(1) = ii(1) + 1 + ell_1
104  if (buf(ind(1), ind(2)) .ne. i) then
105  s = s + 1
106  buf(ind(1), ind(2)) = i
107  end if
108  end do
109  end do
110  end do
111  pre_compute_n(i) = s
112  nb = nb + s
113  end do
114  max_nb = max(max_nb, nb)
115 
116  size_pre_compute = max_nb
117 
118  sll_deallocate_array(buf, error)
119 
120  end subroutine pre_precompute_double_gyroaverage_coeff_polar_splines
121 
122  subroutine precompute_double_gyroaverage_coeff_polar_splines( &
123  r_min, &
124  r_max, &
125  num_cells_r, &
126  num_cells_theta, &
127  rho, &
128  points, &
129  N_points, &
130  size_pre_compute, &
131  pre_compute_index, &
132  pre_compute_coeff_spl)
133  sll_real64, intent(in) :: r_min
134  sll_real64, intent(in) :: r_max
135  sll_int32, intent(in) :: num_cells_r
136  sll_int32, intent(in) :: num_cells_theta
137  sll_real64, dimension(:), intent(in) :: rho
138  sll_real64, dimension(:, :), intent(in) :: points
139  sll_int32, intent(in) :: n_points
140  sll_int32, intent(in) :: size_pre_compute
141  sll_int32, dimension(:, :), intent(out) :: pre_compute_index
142  sll_real, dimension(:), intent(out) :: pre_compute_coeff_spl
143  sll_int32, dimension(:, :), allocatable :: buf
144  sll_int32 ::i, j, k, ell_1, ell_2, ii(2), s, nb, ind(2)
145  sll_real64::val(-1:2, -1:2), eta_star(2), eta(2), delta_eta(2), x(2)
146  sll_int32 ::error
147  sll_real64 :: eta_min(2), eta_max(2)
148  sll_int32 :: nc(2)
149 
150  delta_eta(1) = (r_max - r_min)/real(num_cells_r, f64)
151  delta_eta(2) = 2._f64*sll_p_pi/real(num_cells_theta, f64)
152 
153  nc(1) = num_cells_r
154  nc(2) = num_cells_theta
155 
156  eta_min(1) = r_min
157  eta_min(2) = 0._f64
158  eta_max(1) = r_max
159  eta_max(2) = 2._f64*sll_p_pi
160 
161  sll_allocate(buf(0:num_cells_r + 2, 0:num_cells_theta - 1), error)
162 
163  eta(2) = 0._f64
164 
165  buf = 0
166  nb = 0
167  s = 0
168  val = 0._f64
169  eta(2) = 0._f64
170  do i = 1, nc(1) + 1
171  eta(1) = eta_min(1) + real(i - 1, f64)*delta_eta(1)
172  do k = 1, n_points
173  x(1) = eta(1)*cos(eta(2)) + rho(i)*points(1, k)
174  x(2) = eta(1)*sin(eta(2)) + rho(i)*points(2, k)
175  call sll_s_localize_polar(x, eta_min, eta_max, ii, eta_star, nc)
176 
177  call sll_s_contribution_spl(eta_star, val)
178 
179  val = val*points(3, k)
180 
181  do ell_2 = -1, 2
182  ind(2) = modulo(ii(2) + ell_2, nc(2))
183  do ell_1 = -1, 2
184  ind(1) = ii(1) + 1 + ell_1
185  j = buf(ind(1), ind(2))
186  if (j <= nb) then
187  s = s + 1
188  buf(ind(1), ind(2)) = s
189  pre_compute_coeff_spl(s) = val(ell_1, ell_2)
190  pre_compute_index(1, s) = ind(1)
191  pre_compute_index(2, s) = ind(2)
192  else
193  pre_compute_coeff_spl(j) = pre_compute_coeff_spl(j) + val(ell_1, ell_2)
194  end if
195  end do
196  end do
197  end do
198  nb = s
199  end do
200 
201  sll_deallocate_array(buf, error)
202  return
203  print *, size_pre_compute
204 
205  end subroutine precompute_double_gyroaverage_coeff_polar_splines
206 
207  subroutine compute_contribution_matrix( &
208  num_cells_r, &
209  num_cells_theta, &
210  pre_compute_N, &
211  pre_compute_index, &
212  pre_compute_coeff_spl, &
213  size_pre_compute, &
214  pointer_mat_contribution_circ)
215  sll_int32, intent(in) :: num_cells_r
216  sll_int32, intent(in) :: num_cells_theta
217  sll_int32, dimension(:), intent(in) :: pre_compute_n
218  sll_int32, dimension(:, :), intent(in) :: pre_compute_index
219  sll_real64, dimension(:), intent(in) :: pre_compute_coeff_spl
220  sll_int32, intent(in) :: size_pre_compute
221  sll_real64, dimension(:, :, :), intent(out) :: pointer_mat_contribution_circ
222  sll_int32 :: i
223  sll_int32 :: s
224  sll_int32 :: k
225  sll_int32 :: ii(2)
226  pointer_mat_contribution_circ = 0._f64
227  s = 0
228  do i = 1, num_cells_r + 1
229  do k = 1, pre_compute_n(i)
230  s = s + 1
231  ii(1) = pre_compute_index(1, s) + 1
232  ii(2) = modulo(pre_compute_index(2, s), num_cells_theta) + 1
233  pointer_mat_contribution_circ(ii(2), i, ii(1)) = &
234  pointer_mat_contribution_circ(ii(2), i, ii(1)) + &
235  pre_compute_coeff_spl(s)
236  end do
237  end do
238  return
239  print *, size_pre_compute
240  end subroutine compute_contribution_matrix
241 
242  subroutine compute_d_spl2d( &
243  num_cells_r, &
244  num_cells_theta, &
245  D_spl2D)
246  sll_int32, intent(in) :: num_cells_r
247  sll_int32, intent(in) :: num_cells_theta
248  sll_int32 :: ierr
249  sll_int32 :: nr
250  sll_int32 :: ntheta
251  sll_comp64, dimension(:, :, :), intent(out) :: d_spl2d
252  sll_real64, dimension(:), allocatable, target::dnat, lnat, dper, lper, mper
253  sll_real64, dimension(:), pointer::pointer_dnat, pointer_lnat
254  sll_real64, dimension(:), pointer::pointer_dper, pointer_lper, pointer_mper
255  sll_real64, dimension(:, :), allocatable, target::mat_nat, mat_per
256  sll_real64, dimension(:, :), pointer::pointer_mat_nat, pointer_mat_per
257  sll_real64, dimension(:, :, :), allocatable, target::mat_spl2d_circ
258  sll_real64, dimension(:, :, :), pointer::pointer_mat_spl2d_circ
259  sll_int32 :: j
260  !sll_int32 :: m
261  !sll_comp64 :: exp_comp
262  !sll_real64 :: mode
263  sll_comp64, dimension(:), allocatable :: fft_array
264  sll_real64, dimension(:), allocatable :: buf_fft
265  sll_int32 :: k
266 
267  nr = num_cells_r
268  ntheta = num_cells_theta
269 
270  sll_allocate(dnat(0:nr + 2), ierr)
271  sll_allocate(lnat(0:nr + 2), ierr)
272  sll_allocate(dper(0:ntheta - 1), ierr)
273  sll_allocate(lper(0:ntheta - 1), ierr)
274  sll_allocate(mper(0:ntheta - 1), ierr)
275  sll_allocate(mat_nat(0:nr + 2, 0:nr), ierr)
276  sll_allocate(mat_per(0:ntheta - 1, 0:ntheta - 1), ierr)
277  sll_allocate(mat_spl2d_circ(0:ntheta - 1, 0:nr + 2, 0:nr), ierr)
278  !SLL_ALLOCATE(D_spl2D(0:Ntheta-1,0:Nr+2,0:Nr),error)
279  sll_allocate(buf_fft(4*ntheta + 15), ierr)
280  sll_allocate(fft_array(ntheta), ierr)
281 
282  pointer_dnat => dnat
283  pointer_lnat => lnat
284  pointer_dper => dper
285  pointer_lper => lper
286  pointer_mper => mper
287  pointer_mat_nat => mat_nat
288  pointer_mat_per => mat_per
289  pointer_mat_spl2d_circ => mat_spl2d_circ
290 
291  call sll_s_splcoefnat1d0old(dnat, lnat, nr)
292  call sll_s_splcoefper1d0old(dper, lper, mper, ntheta)
293 
294  call sll_s_compute_splines_coefs_matrix_nat_1d(pointer_mat_nat, pointer_dnat, pointer_lnat, nr)
295  call sll_s_compute_splines_coefs_matrix_per_1d(pointer_mat_per, pointer_dper, pointer_lper, pointer_mper, ntheta)
296  do j = 0, ntheta - 1
297  pointer_mat_spl2d_circ(j, :, :) = pointer_mat_per(0, j)*pointer_mat_nat
298  end do
299 
300  call zffti(ntheta, buf_fft)
301 
302  do k = 0, nr
303  do j = 0, nr + 2
304  fft_array(1:ntheta) = pointer_mat_spl2d_circ(0:ntheta - 1, j, k)*(1._f64, 0._f64)
305  call zfftf(ntheta, fft_array(1:ntheta), buf_fft)
306  d_spl2d(1:ntheta, j + 1, k + 1) = fft_array(1:ntheta)
307  end do
308  end do
309 
310 ! D_spl2D = 0._f64
311 ! do m=0,Ntheta-1
312 ! do j=0,Ntheta-1
313 ! mode=real(-2._f64*sll_p_pi*real(j,f64)*real(m,f64)/real(Ntheta,f64),f64)
314 ! exp_comp = cmplx( cos(mode), sin(mode), kind=f64 )
315 ! D_spl2D(m,:,:) = D_spl2D(m,:,:) + pointer_mat_spl2D_circ(j,:,:)*exp_comp
316 ! enddo
317 ! enddo
318 
319  end subroutine compute_d_spl2d
320 
321  subroutine compute_d_contr( &
322  num_cells_r, &
323  num_cells_theta, &
324  pointer_mat_contribution_circ, &
325  D_contr)
326  sll_int32, intent(in) :: num_cells_r
327  sll_int32, intent(in) :: num_cells_theta
328  sll_real64, dimension(:, :, :), intent(in) :: pointer_mat_contribution_circ
329  sll_int32 :: ierr
330  sll_int32 :: nr
331  sll_int32 :: ntheta
332  sll_comp64, dimension(:, :, :), intent(out) :: d_contr
333  !sll_int32 :: m
334  sll_int32 :: j
335  !sll_real64 :: mode
336  !sll_comp64 :: exp_comp
337  sll_comp64, dimension(:), allocatable :: fft_array
338  sll_real64, dimension(:), allocatable :: buf_fft
339  sll_int32 :: k
340 
341  nr = num_cells_r
342  ntheta = num_cells_theta
343 
344  sll_allocate(buf_fft(4*ntheta + 15), ierr)
345  sll_allocate(fft_array(ntheta), ierr)
346 
347  call zffti(ntheta, buf_fft)
348 
349  do k = 0, nr + 2
350  do j = 0, nr
351  fft_array(1:ntheta) = pointer_mat_contribution_circ(1:ntheta, j + 1, k + 1)*(1._f64, 0._f64)
352  call zfftf(ntheta, fft_array(1:ntheta), buf_fft)
353  d_contr(1:ntheta, j + 1, k + 1) = fft_array(1:ntheta)
354  end do
355  end do
356 
357 ! D_contr = 0._f64
358 ! do m=0,Ntheta-1
359 ! do j=0,Ntheta-1
360 ! mode=real(-2._f64*sll_p_pi*real(j,f64)*real(m,f64)/real(Ntheta,f64),f64)
361 ! exp_comp = cmplx( cos(mode), sin(mode), kind=f64 )
362 ! D_contr(m,:,:) = D_contr(m,:,:) + pointer_mat_contribution_circ(j,:,:)*exp_comp
363 ! enddo
364 ! enddo
365 
366  end subroutine compute_d_contr
367 
368  subroutine compute_double_gyroaverage_matrix( &
369  D_spl2D, &
370  D_contr, &
371  num_cells_r, &
372  num_cells_theta, &
373  mat)
374  sll_comp64, dimension(:, :, :), intent(in) :: d_spl2d
375  sll_comp64, dimension(:, :, :), intent(in) :: d_contr
376  sll_int32, intent(in) :: num_cells_r
377  sll_int32, intent(in) :: num_cells_theta
378  sll_comp64, dimension(:, :, :), intent(out) :: mat
379  sll_int32 :: nr
380  sll_int32 :: ntheta
381  sll_int32 :: m
382  sll_int32 :: i
383  sll_comp64, dimension(:, :), allocatable :: mat_stock1
384  sll_comp64, dimension(:, :), allocatable :: mat_stock2
385  sll_comp64, dimension(:, :), allocatable :: mat_stock3
386  sll_int32 :: ierr
387 
388  nr = num_cells_r
389  ntheta = num_cells_theta
390 
391  sll_allocate(mat_stock1(nr + 1, nr + 1), ierr)
392  sll_allocate(mat_stock2(nr + 1, nr + 1), ierr)
393  sll_allocate(mat_stock3(nr + 1, nr + 1), ierr)
394 
395  mat = (0._f64, 0._f64)
396  do m = 1, ntheta
397  mat_stock1 = (0._f64, 0._f64)
398  mat_stock2 = (0._f64, 0._f64)
399  call sll_s_matrix_product_compf( &
400  d_contr(m, :, :), &
401  nr + 1, &
402  nr + 3, &
403  d_spl2d(m, :, :), &
404  nr + 3, &
405  nr + 1, &
406  mat_stock1)
407  ! call matrix_product_comp( &
408 ! mat_stock1(:,:), &
409 ! Nr+1, &
410 ! Nr+1, &
411 ! mat_stock1(:,:), &
412 ! Nr+1, &
413 ! Nr+1, &
414 ! mat_stock2)
415 !
416  !print *,'mat_stock1 l1',m,mat_stock1(1,1),mat_stock1(1,2)
417  !print *,'mat_stock1 l1',m,mat_stock1(2,1),mat_stock1(2,2)
418 
419  !mat_stock3 = transpose(mat_stock1)
420  mat_stock3 = mat_stock1
421  call sll_s_matrix_product_compf( &
422  mat_stock3, &
423  nr + 1, &
424  nr + 1, &
425  mat_stock1, &
426  nr + 1, &
427  nr + 1, &
428  mat_stock2)
429  ! mat_stock2 = mat_stock1
430 
431 ! if(m==1)then
432 
433  do i = 1, nr + 1
434  mat_stock2(i, i) = mat_stock2(i, i) - (1._f64, 0._f64)
435  end do
436 
437 ! endif
438  mat(m, :, :) = -mat_stock2
439  !mat(m,:,:) = mat_stock2
440 
441  !print *,'mat_stock2 l1',m,mat_stock2(1,1),mat_stock2(1,2)
442  !print *,'mat_stock2 l1',m,mat_stock2(2,1),mat_stock2(2,2)
443 
444  end do
445 
446  end subroutine compute_double_gyroaverage_matrix
447 
448  subroutine compute_double_gyroaverage_matrix_blas( &
449  D_spl2D, &
450  D_contr, &
451  num_cells_r, &
452  num_cells_theta, &
453  mat)
454  sll_comp64, dimension(:, :, :), intent(in) :: d_spl2d
455  sll_comp64, dimension(:, :, :), intent(in) :: d_contr
456  sll_int32, intent(in) :: num_cells_r
457  sll_int32, intent(in) :: num_cells_theta
458  sll_comp64, dimension(:, :, :), intent(out) :: mat
459  sll_int32 :: nr
460  sll_int32 :: ntheta
461  sll_int32 :: m
462  sll_int32 :: i
463  sll_comp64, dimension(:, :), allocatable :: mat_stock1
464  sll_comp64, dimension(:, :), allocatable :: mat_stock2
465 
466  sll_comp64, dimension(:, :), allocatable :: mat_a
467  sll_comp64, dimension(:, :), allocatable :: mat_b
468  sll_comp64, dimension(:, :), allocatable :: mat_c
469 
470  sll_int32 :: ierr
471 
472  nr = num_cells_r
473  ntheta = num_cells_theta
474 
475  sll_allocate(mat_stock1(nr + 1, nr + 1), ierr)
476  sll_allocate(mat_stock2(nr + 1, nr + 1), ierr)
477  sll_allocate(mat_a(nr + 1, nr + 3), ierr)
478  sll_allocate(mat_b(nr + 3, nr + 1), ierr)
479  sll_allocate(mat_c(nr + 1, nr + 1), ierr)
480 
481  mat = (0._f64, 0._f64)
482  do m = 1, ntheta
483  mat_stock1 = (0._f64, 0._f64)
484  mat_a(1:nr + 1, 1:nr + 3) = d_contr(m, :, :)
485  mat_b(1:nr + 3, 1:nr + 1) = d_spl2d(m, :, :)
486  !A MxK
487  !B KxN
488  !C MxN
489  !M = Nr+1, K = Nr+3, N = Nr+1
490  !ALPHA=1 BETA=0
491  !ZGEMM('N','N',M,N,K,ALPHA,A,M,B,K,BETA,C,M)
492  CALL zgemm( &
493  'N', &
494  'N', &
495  nr + 1, &
496  nr + 1, &
497  nr + 3, &
498  1._f64, &
499  mat_a, &
500  nr + 1, &
501  mat_b, &
502  nr + 3, &
503  0._f64, &
504  mat_c, &
505  nr + 1)
506 
507 ! call matrix_product_comp( &
508 ! D_contr(m,:,:), &
509 ! Nr+1, &
510 ! Nr+3, &
511 ! D_spl2D(m,:,:), &
512 ! Nr+3, &
513 ! Nr+1, &
514 ! mat_stock1)
515 
516  mat_stock2 = (0._f64, 0._f64)
517  mat_stock1(:, :) = mat_c(1:nr + 1, 1:nr + 1)
518  CALL zgemm( &
519  'T', &
520  'N', &
521  nr + 1, &
522  nr + 1, &
523  nr + 1, &
524  1._f64, &
525  mat_c, &
526  nr + 1, &
527  mat_stock1(:, :), &
528  nr + 1, &
529  0._f64, &
530  mat_stock2, &
531  nr + 1)
532 
533 ! call matrix_product_comp( &
534 ! transpose(mat_stock1(:,:)), &
535 ! Nr+1, &
536 ! Nr+1, &
537 ! mat_stock1(:,:), &
538 ! Nr+1, &
539 ! Nr+1, &
540 ! mat_stock2)
541  do i = 1, nr + 1
542  mat_stock2(i, i) = mat_stock2(i, i) - (1._f64, 0._f64)
543  end do
544  mat(m, :, :) = -mat_stock2
545  end do
546 
547  end subroutine compute_double_gyroaverage_matrix_blas
548 
549  subroutine compute_double_gyroaverage_matrix_blas_real( &
550  D_spl2D, &
551  D_contr, &
552  num_cells_r, &
553  num_cells_theta, &
554  mat)
555  sll_real64, dimension(:, :, :), intent(in) :: d_spl2d
556  sll_real64, dimension(:, :, :), intent(in) :: d_contr
557  sll_int32, intent(in) :: num_cells_r
558  sll_int32, intent(in) :: num_cells_theta
559  sll_real64, dimension(:, :, :), intent(out) :: mat
560  sll_int32 :: nr
561  sll_int32 :: ntheta
562  sll_int32 :: m
563  sll_int32 :: i
564  sll_real64, dimension(:, :), allocatable :: mat_stock1
565  sll_real64, dimension(:, :), allocatable :: mat_stock2
566 
567  sll_real64, dimension(:, :), allocatable :: mat_a
568  sll_real64, dimension(:, :), allocatable :: mat_b
569  sll_real64, dimension(:, :), allocatable :: mat_c
570 
571  sll_int32 :: ierr
572 
573  nr = num_cells_r
574  ntheta = num_cells_theta
575 
576  sll_allocate(mat_stock1(nr + 1, nr + 1), ierr)
577  sll_allocate(mat_stock2(nr + 1, nr + 1), ierr)
578  sll_allocate(mat_a(nr + 1, nr + 3), ierr)
579  sll_allocate(mat_b(nr + 3, nr + 1), ierr)
580  sll_allocate(mat_c(nr + 1, nr + 1), ierr)
581 
582  mat = 0._f64
583  do m = 1, ntheta
584  mat_stock1 = 0._f64
585  mat_a(1:nr + 1, 1:nr + 3) = d_contr(m, :, :)
586  mat_b(1:nr + 3, 1:nr + 1) = d_spl2d(m, :, :)
587  !A MxK
588  !B KxN
589  !C MxN
590  !M = Nr+1, K = Nr+3, N = Nr+1
591  !ALPHA=1 BETA=0
592  !ZGEMM('N','N',M,N,K,ALPHA,A,M,B,K,BETA,C,M)
593  CALL dgemm( &
594  'N', &
595  'N', &
596  nr + 1, &
597  nr + 1, &
598  nr + 3, &
599  1._f64, &
600  mat_a, &
601  nr + 1, &
602  mat_b, &
603  nr + 3, &
604  0._f64, &
605  mat_c, &
606  nr + 1)
607 
608 ! call matrix_product_comp( &
609 ! D_contr(m,:,:), &
610 ! Nr+1, &
611 ! Nr+3, &
612 ! D_spl2D(m,:,:), &
613 ! Nr+3, &
614 ! Nr+1, &
615 ! mat_stock1)
616 
617  mat_stock2 = 0._f64
618  mat_stock1(:, :) = mat_c(1:nr + 1, 1:nr + 1)
619  CALL dgemm( &
620  'T', &
621  'N', &
622  nr + 1, &
623  nr + 1, &
624  nr + 1, &
625  1._f64, &
626  mat_c, &
627  nr + 1, &
628  mat_stock1(:, :), &
629  nr + 1, &
630  0._f64, &
631  mat_stock2, &
632  nr + 1)
633 
634 ! call matrix_product_comp( &
635 ! transpose(mat_stock1(:,:)), &
636 ! Nr+1, &
637 ! Nr+1, &
638 ! mat_stock1(:,:), &
639 ! Nr+1, &
640 ! Nr+1, &
641 ! mat_stock2)
642  do i = 1, nr + 1
643  mat_stock2(i, i) = mat_stock2(i, i) - 1._f64
644  end do
645  mat(m, :, :) = -mat_stock2
646  end do
647 
648  end subroutine compute_double_gyroaverage_matrix_blas_real
649 
650  subroutine sll_s_compute_qns_matrix_polar_splines( &
651  r_min, &
652  r_max, &
653  num_cells_r, &
654  num_cells_theta, &
655  rho, &
656  points, &
657  N_points, &
658  mat)
659  sll_real64, intent(in) :: r_min
660  sll_real64, intent(in) :: r_max
661  sll_int32, intent(in) :: num_cells_r
662  sll_int32, intent(in) :: num_cells_theta
663  sll_real64, dimension(:), intent(in) :: rho
664  sll_real64, dimension(:, :), intent(in) :: points
665  sll_int32, intent(in) :: n_points
666  sll_comp64, dimension(:, :, :), intent(out) :: mat
667  sll_int32, dimension(:), allocatable :: pre_compute_n
668  sll_int32 :: size_pre_compute
669  sll_int32, dimension(:, :), allocatable :: pre_compute_index
670  sll_real64, dimension(:), allocatable :: pre_compute_coeff_spl
671  sll_real64, dimension(:, :, :), allocatable :: pointer_mat_contribution_circ
672  sll_comp64, dimension(:, :, :), allocatable :: d_spl2d
673  sll_comp64, dimension(:, :, :), allocatable :: d_contr
674  sll_int32 :: ierr
675  !sll_real64, dimension(:,:,:), allocatable :: real_mat
676  !sll_real64, dimension(:,:,:), allocatable :: real_D_spl2D
677  !sll_real64, dimension(:,:,:), allocatable :: real_D_contr
678 
679  sll_allocate(pre_compute_n(num_cells_r + 1), ierr)
680 
681  call pre_precompute_double_gyroaverage_coeff_polar_splines( &
682  r_min, &
683  r_max, &
684  num_cells_r, &
685  num_cells_theta, &
686  rho, &
687  points, &
688  n_points, &
689  pre_compute_n, &
690  size_pre_compute)
691 
692  !print *,'#size_pre_compute=',size_pre_compute
693 
694  sll_allocate(pre_compute_index(2, size_pre_compute), ierr)
695  sll_allocate(pre_compute_coeff_spl(size_pre_compute), ierr)
696 
697  call precompute_double_gyroaverage_coeff_polar_splines( &
698  r_min, &
699  r_max, &
700  num_cells_r, &
701  num_cells_theta, &
702  rho, &
703  points, &
704  n_points, &
705  size_pre_compute, &
706  pre_compute_index, &
707  pre_compute_coeff_spl)
708 
709  !print *,'#precompute_double_gyroaverage_coeff_polar_splines done'
710 
711  sll_allocate(pointer_mat_contribution_circ(num_cells_theta, num_cells_r + 1, num_cells_r + 3), ierr)
712 
713  call compute_contribution_matrix( &
714  num_cells_r, &
715  num_cells_theta, &
716  pre_compute_n, &
717  pre_compute_index, &
718  pre_compute_coeff_spl, &
719  size_pre_compute, &
720  pointer_mat_contribution_circ)
721 
722  !print *,'#compute_contribution_matrix done'
723 
724  sll_allocate(d_spl2d(num_cells_theta, num_cells_r + 3, num_cells_r + 1), ierr)
725 
726  call compute_d_spl2d( &
727  num_cells_r, &
728  num_cells_theta, &
729  d_spl2d)
730 
731  !print *,'#compute_D_spl2D done'
732 
733  sll_allocate(d_contr(num_cells_theta, num_cells_r + 1, num_cells_r + 3), ierr)
734 
735  call compute_d_contr( &
736  num_cells_r, &
737  num_cells_theta, &
738  pointer_mat_contribution_circ, &
739  d_contr)
740 
741  !print *,'#compute_D_contr done'
742 
743  !print *,'#complex values of D_spl2D',minval(aimag(D_spl2D)),maxval(aimag(D_spl2D))
744  !print *,'#complex values of D_contr',minval(aimag(D_contr)),maxval(aimag(D_contr))
745  !print *,'#real values of D_spl2D',minval(real(D_spl2D)),maxval(real(D_spl2D))
746  !print *,'#real values of D_contr',minval(real(D_contr)),maxval(real(D_contr))
747 
748  !print *,'#begin compute_double_gyroaverage_matrix'
749 
750  call compute_double_gyroaverage_matrix( &
751  d_spl2d, &
752  d_contr, &
753  num_cells_r, &
754  num_cells_theta, &
755  mat)
756  !print *,'#end compute_double_gyroaverage_matrix'
757 
758 ! call compute_double_gyroaverage_matrix_blas( &
759 ! D_spl2D, &
760 ! D_contr, &
761 ! num_cells_r, &
762 ! num_cells_theta, &
763 ! mat)
764 
765 ! SLL_ALLOCATE(real_D_spl2D(num_cells_theta,num_cells_r+3,num_cells_r+1), ierr)
766 ! SLL_ALLOCATE(real_D_contr(num_cells_theta,num_cells_r+1,num_cells_r+3), ierr)
767 ! SLL_ALLOCATE(real_mat(num_cells_theta,num_cells_r+1,num_cells_r+1), ierr)
768 !
769 ! real_D_spl2D = real(D_spl2D,f64)
770 ! real_D_contr = real(D_contr,f64)
771 !
772 ! real_mat = 0._f64
773 !
774 !
775 !
776 ! call compute_double_gyroaverage_matrix_blas_real( &
777 ! real_D_spl2D, &
778 ! real_D_contr, &
779 ! num_cells_r, &
780 ! num_cells_theta, &
781 ! real_mat)
782 !
783 ! mat = real_mat*(1._f64,0._f64)
784 
785  !print *,'#compute_double_gyroaverage_matrix'
786 
787  end subroutine sll_s_compute_qns_matrix_polar_splines
788 
789  subroutine sll_s_compute_qns_inverse_polar_splines( &
790  mat, &
791  lambda, &
792  num_cells_r, &
793  num_cells_theta)
794  sll_comp64, dimension(:, :, :), intent(inout) :: mat
795  sll_real64, dimension(:), intent(in) :: lambda
796  sll_int32, intent(in) :: num_cells_r
797  sll_int32, intent(in) :: num_cells_theta
798  sll_int32 :: ierr
799  sll_int32, dimension(:), allocatable :: ipiv
800  sll_comp64, dimension(:), allocatable :: work
801  sll_int32 :: info
802  sll_int32 :: i
803  sll_int32 :: m
804 
805  sll_allocate(ipiv(num_cells_r + 1), ierr)
806  sll_allocate(work((num_cells_r + 1)**2), ierr)
807  do i = 1, num_cells_r + 1
808  do m = 1, num_cells_theta
809  mat(m, i, i) = mat(m, i, i) + lambda(i)
810  end do
811  end do
812 
813  do m = 1, num_cells_theta
814  call zgetrf( &
815  num_cells_r + 1, &
816  num_cells_r + 1, &
817  mat(m, :, :), &
818  num_cells_r + 1, &
819  ipiv, &
820  info)
821  call zgetri( &
822  num_cells_r + 1, &
823  mat(m, :, :), &
824  num_cells_r + 1, &
825  ipiv, &
826  work, &
827  (num_cells_r + 1)**2, &
828  info)
829  end do
830 
831  print *, '#here is INFO'
832  print *, info
833 
834  end subroutine sll_s_compute_qns_inverse_polar_splines
835 
836  subroutine sll_s_solve_qns_polar_splines( &
837  mat, &
838  phi, &
839  num_cells_r, &
840  num_cells_theta)
841  sll_comp64, dimension(:, :, :), intent(in) :: mat
842  sll_real64, dimension(:, :), intent(inout) :: phi
843  sll_comp64, dimension(:, :), allocatable :: phi_comp
844  sll_comp64, dimension(:, :), allocatable :: phi_old
845  sll_real64, dimension(:), allocatable::buf_fft
846  sll_int32, intent(in) :: num_cells_r
847  sll_int32, intent(in) :: num_cells_theta
848  sll_int32 :: nr
849  sll_int32 :: ntheta
850  sll_int32 :: m
851  sll_int32 :: i
852  sll_int32 :: j
853  sll_int32 :: ierr
854  sll_comp64 :: result
855 
856  nr = num_cells_r
857  ntheta = num_cells_theta
858 
859  sll_allocate(phi_comp(1:nr + 1, 1:ntheta), ierr)
860  sll_allocate(phi_old(1:nr + 1, 1:ntheta), ierr)
861  sll_allocate(buf_fft(1:4*ntheta + 15), ierr)
862 
863  ! FFT(PHI)
864  phi_comp = phi*(1._f64, 0._f64)
865  call zffti(ntheta, buf_fft)
866  do i = 1, nr + 1
867  call zfftf(ntheta, phi_comp(i, :), buf_fft)
868  end do
869 
870  ! Produit matrice/vecteur
871  phi_old = phi_comp
872  do m = 1, ntheta
873  do i = 1, nr + 1
874  result = (0._f64, 0._f64)
875  do j = 1, nr + 1
876  result = result + mat(m, i, j)*phi_old(j, m)
877  end do
878  phi_comp(i, m) = result
879  end do
880  end do
881 
882  ! FFT^-1
883  do i = 1, nr + 1
884  call zfftb(ntheta, phi_comp(i, :), buf_fft)
885  end do
886  phi = real(phi_comp, f64)/real(ntheta, f64)
887 
888  ! Sorties
889  ! print *,"phi_min : ",minval(phi)
890  ! print *,"phi_max : ",maxval(phi)
891 
892  ! Deallocate
893  sll_deallocate_array(phi_comp, ierr)
894  sll_deallocate_array(phi_old, ierr)
895  sll_deallocate_array(buf_fft, ierr)
896 
897  end subroutine sll_s_solve_qns_polar_splines
898 
899 !PN subroutine precompute_matrix( &
900 !PN Ti, &
901 !PN r_min, &
902 !PN r_max, &
903 !PN num_cells_r, &
904 !PN num_cells_theta, &
905 !PN mu_points, &
906 !PN mu_weights, &
907 !PN N_mu, &
908 !PN N_points, &
909 !PN mat)
910 !PN sll_real64,dimension(:),intent(in) :: Ti
911 !PN sll_real64, intent(in) :: r_min
912 !PN sll_real64, intent(in) :: r_max
913 !PN sll_int32, intent(in) :: num_cells_r
914 !PN sll_int32, intent(in) :: num_cells_theta
915 !PN sll_real64, dimension(:), intent(in) :: mu_points
916 !PN sll_real64, dimension(:), intent(in) :: mu_weights
917 !PN sll_int32, intent(in) :: N_mu
918 !PN sll_int32, intent(in) :: N_points
919 !PN sll_real64, dimension(:,:,:), intent(out) :: mat
920 !PN
921 !PN
922 !PN end subroutine precompute_matrix
923 
924 end module sll_m_qn_2d_polar_precompute
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_qn_2d_polar_precompute
Definition: sll_m_qn_2d_polar_precompute.F90:1
sll_m_qn_2d_polar_precompute::sll_s_compute_qns_matrix_polar_splines
subroutine, public sll_s_compute_qns_matrix_polar_splines(r_min, r_max, num_cells_r, num_cells_theta, rho, points, N_points, mat)
Definition: sll_m_qn_2d_polar_precompute.F90:659
sll_m_qn_2d_polar_precompute::compute_d_contr
subroutine compute_d_contr(num_cells_r, num_cells_theta, pointer_mat_contribution_circ, D_contr)
Definition: sll_m_qn_2d_polar_precompute.F90:326
sll_m_qn_2d_polar_precompute::precompute_double_gyroaverage_coeff_polar_splines
subroutine precompute_double_gyroaverage_coeff_polar_splines(r_min, r_max, num_cells_r, num_cells_theta, rho, points, N_points, size_pre_compute, pre_compute_index, pre_compute_coeff_spl)
Definition: sll_m_qn_2d_polar_precompute.F90:133
sll_m_qn_2d_polar_precompute::compute_d_spl2d
subroutine compute_d_spl2d(num_cells_r, num_cells_theta, D_spl2D)
Definition: sll_m_qn_2d_polar_precompute.F90:246
sll_m_qn_2d_polar_precompute::sll_s_solve_qns_polar_splines
subroutine, public sll_s_solve_qns_polar_splines(mat, phi, num_cells_r, num_cells_theta)
Definition: sll_m_qn_2d_polar_precompute.F90:841
sll_m_qn_2d_polar_precompute::compute_contribution_matrix
subroutine compute_contribution_matrix(num_cells_r, num_cells_theta, pre_compute_N, pre_compute_index, pre_compute_coeff_spl, size_pre_compute, pointer_mat_contribution_circ)
Definition: sll_m_qn_2d_polar_precompute.F90:215
sll_m_qn_2d_polar_precompute::compute_double_gyroaverage_matrix_blas
subroutine compute_double_gyroaverage_matrix_blas(D_spl2D, D_contr, num_cells_r, num_cells_theta, mat)
Definition: sll_m_qn_2d_polar_precompute.F90:454
sll_m_qn_2d_polar_precompute::compute_double_gyroaverage_matrix
subroutine compute_double_gyroaverage_matrix(D_spl2D, D_contr, num_cells_r, num_cells_theta, mat)
Definition: sll_m_qn_2d_polar_precompute.F90:374
sll_m_qn_2d_polar_precompute::sll_s_compute_qns_inverse_polar_splines
subroutine, public sll_s_compute_qns_inverse_polar_splines(mat, lambda, num_cells_r, num_cells_theta)
Definition: sll_m_qn_2d_polar_precompute.F90:794
sll_m_qn_2d_polar_precompute::compute_double_gyroaverage_matrix_blas_real
subroutine compute_double_gyroaverage_matrix_blas_real(D_spl2D, D_contr, num_cells_r, num_cells_theta, mat)
Definition: sll_m_qn_2d_polar_precompute.F90:555
sll_m_qn_2d_polar_precompute::pre_precompute_double_gyroaverage_coeff_polar_splines
subroutine pre_precompute_double_gyroaverage_coeff_polar_splines(r_min, r_max, num_cells_r, num_cells_theta, rho, points, N_points, pre_compute_N, size_pre_compute)
Definition: sll_m_qn_2d_polar_precompute.F90:53
sll_m_qn_2d_polar
Definition: sll_m_qn_2d_polar.F90:18
sll_m_qn_2d_polar::sll_s_compute_splines_coefs_matrix_per_1d
subroutine, public sll_s_compute_splines_coefs_matrix_per_1d(mat, dper, lper, mper, N)
Definition: sll_m_qn_2d_polar.F90:169
sll_m_qn_2d_polar::sll_s_contribution_spl
subroutine, public sll_s_contribution_spl(x, val)
Definition: sll_m_qn_2d_polar.F90:1283
sll_m_qn_2d_polar::sll_s_matrix_product_compf
subroutine, public sll_s_matrix_product_compf(A, NxA, NyA, B, NxB, NyB, prod)
Definition: sll_m_qn_2d_polar.F90:257
sll_m_qn_2d_polar::sll_s_localize_polar
subroutine, public sll_s_localize_polar(x, eta_min, eta_max, ii, eta, N)
Definition: sll_m_qn_2d_polar.F90:1100
sll_m_qn_2d_polar::sll_s_splcoefper1d0old
subroutine, public sll_s_splcoefper1d0old(dper, lper, mper, N)
Definition: sll_m_qn_2d_polar.F90:1503
sll_m_qn_2d_polar::sll_s_compute_splines_coefs_matrix_nat_1d
subroutine, public sll_s_compute_splines_coefs_matrix_nat_1d(mat, dnat, lnat, N)
Definition: sll_m_qn_2d_polar.F90:124
sll_m_qn_2d_polar::sll_s_splcoefnat1d0old
subroutine, public sll_s_splcoefnat1d0old(dnat, lnat, N)
Definition: sll_m_qn_2d_polar.F90:1361
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:39 for Semi-Lagrangian Library by   doxygen 1.9.1