Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_particle_sampling.F90
Go to the documentation of this file.
1 
5 module sll_m_particle_sampling
6  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7 #include "sll_errors.h"
8 #include "sll_memory.h"
9 #include "sll_assert.h"
10 #include "sll_working_precision.h"
11 
12  use sll_m_collective, only: &
13  sll_o_collective_allreduce, &
14  sll_v_world_collective
15 
16  use sll_m_control_variate, only : &
17  sll_t_control_variate
18 
19  use sll_m_initial_distribution, only: &
20  sll_c_distribution_params, &
21  sll_t_params_cos_gaussian, &
22  sll_t_params_cos_gaussian_screwpinch, &
23  sll_t_params_noise_gaussian
24 
25  use sll_m_particle_group_base, only: &
26  sll_c_particle_group_base
27 
28  use sll_m_prob, only: &
29  sll_s_normal_cdf_inv
30 
31  use sll_m_sobol, only: &
32  sll_s_i8_sobol
33 
34  use sll_m_mapping_3d, only: &
35  sll_t_mapping_3d
36 
37  use mpi, only: &
38  mpi_sum
39 
40  implicit none
41 
42  public :: &
43  sll_t_particle_sampling, &
44  sll_s_particle_sampling_randomized_weights
45 
46  private
47  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
48 
49  ! Internal parameter to distinguish how to draw
50  sll_int32, parameter :: sll_p_random_numbers = 0
51  sll_int32, parameter :: sll_p_sobol_numbers = 1
52 
53  sll_int32, parameter :: sll_p_standard = 0
54  sll_int32, parameter :: sll_p_symmetric_all = 1
55  sll_int32, parameter :: sll_p_symmetric_negative = 2
56  sll_int32, parameter :: sll_p_uniformx_negative = 3
57 
59  type :: sll_t_particle_sampling
60  sll_int32 :: symmetric
61  sll_int32 :: random_numbers
62  logical :: uniform = .false.
63  sll_int32, allocatable :: random_seed(:)
64  sll_int32, allocatable :: random_seed_start(:)
65  sll_int64 :: sobol_seed
66  sll_int64 :: sobol_seed_start
67  logical :: inverse = .false.
68  logical :: xiprofile = .false.
69 
70  logical :: delta_perturb = .false.
71  sll_real64 :: eps(3) = [0.6, 0.6, 0.6]
72  sll_real64 :: a(3) = [0.0, 0.0, 0.0]
73 
74  contains
75  procedure :: init => init_particle_sampling
76  procedure :: sample => sample_particle_sampling
77  procedure :: sample_cv => sample_cv_particle_sampling
78  procedure :: free => free_particle_sampling
79  procedure :: reset_seed_jump
80 
81  end type sll_t_particle_sampling
82 
83 
84 contains
85 
87  subroutine free_particle_sampling( self )
88  class( sll_t_particle_sampling ), intent( inout ) :: self
89 
90  if (allocated(self%random_seed)) deallocate( self%random_seed )
91 
92  end subroutine free_particle_sampling
93 
95  subroutine init_particle_sampling( self, sampling_type, dims, n_particles_local, rank, delta_perturb, delta_eps )
96  class( sll_t_particle_sampling ), intent( out ) :: self
97  character(len=*), intent( in ) :: sampling_type
98  sll_int32, intent( in ) :: dims(:)
99  sll_int32, intent( inout ) :: n_particles_local
100  sll_int32, optional, intent( in ) :: rank
101  logical, optional, intent( in ) :: delta_perturb
102  sll_real64, optional, intent( in ) :: delta_eps(6)
103  ! local variables
104  sll_int32 :: prank
105  sll_int32 :: ncopies
106  sll_int32 :: np, j, rnd_seed_size
107 
108  if( present(delta_perturb) )then
109  self%delta_perturb = delta_perturb
110  if( delta_perturb ) then
111  self%eps = delta_eps(1:3)
112  self%a = delta_eps(4:6)
113  end if
114  end if
115 
116  prank = 0
117  if( present(rank)) prank = rank
118 
119  select case( trim(sampling_type) )
120  case( "particle_sampling_random" )
121  self%symmetric = 0
122  self%random_numbers = sll_p_random_numbers
123  case( "particle_sampling_random_xiprofile" )
124  self%symmetric = 0
125  self%random_numbers = sll_p_random_numbers
126  self%xiprofile = .true.
127  case( "particle_sampling_random_inverse" )
128  self%symmetric = 0
129  self%random_numbers = sll_p_random_numbers
130  self%inverse = .true.
131  case( "particle_sampling_random_inverse_xiprofile" )
132  self%symmetric = 0
133  self%random_numbers = sll_p_random_numbers
134  self%inverse = .true.
135  self%xiprofile = .true.
136  case( "particle_sampling_sobol" )
137  self%symmetric = 0
138  self%random_numbers = sll_p_sobol_numbers
139  case( "particle_sampling_sobol_xiprofile" )
140  self%symmetric = 0
141  self%random_numbers = sll_p_sobol_numbers
142  self%xiprofile = .true.
143  case( "particle_sampling_sobol_inverse" )
144  self%symmetric = 0
145  self%random_numbers = sll_p_sobol_numbers
146  self%inverse = .true.
147  case( "particle_sampling_sobol_inverse_xiprofile" )
148  self%symmetric = 0
149  self%random_numbers = sll_p_sobol_numbers
150  self%inverse = .true.
151  self%xiprofile = .true.
152  case( "particle_sampling_random_symmetric" )
153  self%symmetric = 1
154  self%random_numbers = sll_p_random_numbers
155  case( "particle_sampling_random_symmetric_xiprofile" )
156  self%symmetric = 1
157  self%random_numbers = sll_p_random_numbers
158  self%xiprofile = .true.
159  case( "particle_sampling_random_symmetric_inverse" )
160  self%symmetric = 1
161  self%random_numbers = sll_p_random_numbers
162  self%inverse = .true.
163  case( "particle_sampling_random_symmetric_inverse_xiprofile" )
164  self%symmetric = 1
165  self%random_numbers = sll_p_random_numbers
166  self%inverse = .true.
167  self%xiprofile = .true.
168  case( "particle_sampling_sobol_symmetric" )
169  self%symmetric = 1
170  self%random_numbers = sll_p_sobol_numbers
171  case( "particle_sampling_sobol_symmetric_xiprofile" )
172  self%symmetric = 1
173  self%random_numbers = sll_p_sobol_numbers
174  self%xiprofile = .true.
175  case( "particle_sampling_sobol_symmetric_inverse" )
176  self%symmetric = 1
177  self%random_numbers = sll_p_sobol_numbers
178  self%inverse = .true.
179  case( "particle_sampling_sobol_symmetric_inverse_xiprofile" )
180  self%symmetric = 1
181  self%random_numbers = sll_p_sobol_numbers
182  self%inverse = .true.
183  self%xiprofile = .true.
184  case( "particle_sampling_sobol_symmetric_uniform" )
185  self%symmetric = 1
186  self%uniform = .true.
187  self%random_numbers = sll_p_sobol_numbers
188  case( "particle_sampling_random_negative" )
189  self%symmetric = 2
190  self%random_numbers = sll_p_random_numbers
191  case( "particle_sampling_sobol_negative" )
192  self%symmetric = 2
193  self%random_numbers = sll_p_sobol_numbers
194  case( "particle_sampling_random_unix" )
195  self%symmetric = 3
196  self%random_numbers = sll_p_random_numbers
197  case( "particle_sampling_sobol_unix" )
198  self%symmetric = 3
199  self%random_numbers = sll_p_sobol_numbers
200  case default
201  sll_error("init_particle_sampling","Sampling type not implemented")
202  end select
203 
204  ! Make sure that the particle number is conforming with symmetric sampling (if necessary)
205  if (self%symmetric == 1 .or. self%symmetric == 3 ) then
206  ncopies = 2**(sum(dims))
207  np = modulo(n_particles_local,ncopies)
208  if ( np .ne. 0 ) then
209  n_particles_local = n_particles_local + ncopies - np
210  !print*, n_particles_local, np
211  print*, 'Warning: Number of particles has been increased so that the particle number is compatable with symmetric sampling.'
212  end if
213  elseif (self%symmetric == 2 ) then
214  ncopies = 2
215  np = modulo(n_particles_local,ncopies)
216  if ( np .ne. 0 ) then
217  n_particles_local = n_particles_local + ncopies - np
218  !print*, n_particles_local, np
219  print*, 'Warning: Number of particles has been increased so that the particle number is compatable with symmetric sampling.'
220  end if
221  else
222  ncopies = 1
223  end if
224  ! Set random numbers to default values (overwrite manually if other values are required)
225  select case ( self%random_numbers )
226  case( sll_p_random_numbers )
227  call random_seed(size=rnd_seed_size)
228  allocate(self%random_seed(1:rnd_seed_size))
229  allocate(self%random_seed_start(1:rnd_seed_size))
230  do j=1, rnd_seed_size
231  self%random_seed(j) = (-1)**j*(100 + 15*j)*(2*prank + 1)
232  end do
233  self%random_seed_start = self%random_seed
234  case( sll_p_sobol_numbers )
235  self%sobol_seed = int(10 + prank*n_particles_local/ncopies, 8)
236  self%sobol_seed_start = self%sobol_seed
237  end select
238 
239 
240  end subroutine init_particle_sampling
241 
242  subroutine reset_seed_jump( self, jump )
243  class( sll_t_particle_sampling ), intent( inout ) :: self
244  sll_int32, intent( in ) :: jump
245 
246  ! Set random numbers to default values (overwrite manually if other values are required)
247  select case ( self%random_numbers )
248  case( sll_p_random_numbers )
249  self%random_seed = self%random_seed_start+jump
250  case( sll_p_sobol_numbers )
251  self%sobol_seed = self%sobol_seed_start+int(jump,8)
252  end select
253 
254  end subroutine reset_seed_jump
255 
258  subroutine sample_cv_particle_sampling( self, particle_group, params, xmin, Lx, control_variate, time, map, lindf )
259  class( sll_t_particle_sampling ), intent( inout ) :: self
260  class(sll_c_particle_group_base), target, intent(inout) :: particle_group
261  class( sll_c_distribution_params ), target, intent( inout ) :: params
262  sll_real64, intent(in) :: xmin(:)
263  sll_real64, intent(in) :: lx(:)
264  class(sll_t_control_variate), intent(in) :: control_variate
265  sll_real64, optional, intent(in) :: time
266  type(sll_t_mapping_3d), optional :: map
267  logical, optional :: lindf
268  !local variables
269  sll_real64 :: vi(3), xi(3), x(3), wi(particle_group%n_weights)
270  sll_int32 :: i_part, counter
271  sll_real64 :: time0, rx(3), q, m, g0, df_weight
272  logical :: delta_perturb, lindf0
273 
274  time0 = 0.0_f64
275  lindf0 = .false.
276  if( present(time) ) time0 = time
277  if( present(lindf) ) lindf0 = lindf
278 
279  delta_perturb = self%delta_perturb
280  self%delta_perturb = .false.
281 
282  if( present(map)) then
283  ! First sample the particles and set weights for full f
284  call self%sample( particle_group, params, xmin, lx, map )
285  else
286  ! First sample the particles and set weights for full f
287  call self%sample( particle_group, params, xmin, lx )
288  end if
289 
290  counter = 0
291  ! Fill g0 with value of initial distribution at initial positions (g0)
292  ! and df_weight with the delta f weights
293  do i_part = 1, particle_group%n_particles
294  q = particle_group%species%q
295  m = particle_group%species%m
296  xi = particle_group%get_x( i_part )
297  vi = particle_group%get_v( i_part )
298  wi = particle_group%get_weights( i_part )
299  ! TODO: Distinguish here between different initial sampling distributions
300  if( self%xiprofile ) then
301  if( present(map)) then
302  if( self%inverse ) then
303  g0 = params%eval_v_density( vi(1:params%dims(2)), xi(1:params%dims(1)), m=particle_group%species%m )/map%volume
304  else
305  g0 = params%eval_v_density( vi(1:params%dims(2)), xi(1:params%dims(1)), m=particle_group%species%m )/map%jacobian(xi)
306  end if
307  rx = xi
308  rx(1) = map%get_x1(xi) + m*vi(2)/q
309  rx(1) = (rx(1) - xmin(1))/lx(1)
310  df_weight = control_variate%update_df_weight( rx, vi, time0, wi(1), g0 )
311  else
312  x = (xi - xmin)/lx
313  g0 = params%eval_v_density( vi(1:params%dims(2)), x(1:params%dims(1)), m=particle_group%species%m )/product(lx)
314  rx = x
315  rx(1) = xi(1) + m*vi(2)/q
316  rx(1) = (rx(1) - xmin(1))/lx(1)
317  df_weight = control_variate%update_df_weight( rx, vi, time0, wi(1), g0 )
318  end if
319  else
320  if( present(map)) then
321  x = xi
322  g0 = params%eval_v_density( vi(1:params%dims(2)), xi(1:params%dims(1)), m=particle_group%species%m )/map%jacobian(xi)
323  else
324  x = (xi - xmin)/lx
325  g0 = params%eval_v_density( vi(1:params%dims(2)), xi(1:params%dims(1)), m=particle_group%species%m )/product(lx)
326  end if
327  df_weight = control_variate%update_df_weight( xi, vi, time0, wi(1), g0 )
328  if(delta_perturb) then
329  if( particle_group%species%q < 0._f64) then
330  if( x(1) > 0.5_f64-self%eps(1) .and. x(1) < 0.5_f64+self%eps(1) )then
331  if( x(2) > 0.5_f64-self%eps(2) .and. x(2) < 0.5_f64+self%eps(2) ) then
332  if(xi(3) > 0.5_f64-self%eps(3) .and. x(3) < 0.5_f64+self%eps(3) )then
333  vi = vi + self%a
334  counter = counter + 1
335  end if
336  end if
337  end if
338  end if
339  end if
340  call particle_group%set_v(i_part, vi)
341  end if
342  if(lindf0) then
343  wi(1) = df_weight
344  else
345  wi(2) = g0
346  wi(3) = df_weight
347  end if
348  call particle_group%set_weights( i_part, wi )
349  end do
350  if(delta_perturb) then
351  print*, 'delta_perturb_cv: ', self%eps, self%a, 'shifted particles:', counter
352  end if
353  end subroutine sample_cv_particle_sampling
354 
355 
356 !!$ !> Sample with control variate (we assume that the particle weights are given in the following order:
357 !!$ !> (full f weight, value of initial distribution at time 0, delta f weights)
358 !!$ subroutine sample_lindf_particle_sampling( self, particle_group, params, xmin, Lx, control_variate, time, map )
359 !!$ class( sll_t_particle_sampling ), intent( inout ) :: self !< particle sampling object
360 !!$ class(sll_c_particle_group_base), target, intent(inout) :: particle_group !< particle group
361 !!$ class( sll_c_distribution_params ), target, intent( inout ) :: params !< parameters for initial distribution
362 !!$ sll_real64, intent(in) :: xmin(:) !< lower bound of the domain
363 !!$ sll_real64, intent(in) :: Lx(:) !< length of the domain.
364 !!$ class(sll_t_control_variate), intent(in) :: control_variate !< PIC control variate
365 !!$ sll_real64, optional, intent(in) :: time !< initial time (default: 0)
366 !!$ type(sll_t_mapping_3d), optional :: map !< coordinate transformation
367 !!$ !local variables
368 !!$ sll_real64 :: vi(3), xi(3), x(3), wi(1)
369 !!$ sll_int32 :: i_part, counter
370 !!$ sll_real64 :: time0, Rx(3), q, m
371 !!$ logical :: delta_perturb
372 !!$
373 !!$ time0 = 0.0_f64
374 !!$ if( present(time) ) time0 = time
375 !!$
376 !!$ delta_perturb = self%delta_perturb
377 !!$ self%delta_perturb = .false.
378 !!$
379 !!$ if( present(map)) then
380 !!$ ! First sample the particles and set weights for full f
381 !!$ call self%sample( particle_group, params, xmin, Lx, map )
382 !!$ else
383 !!$ ! First sample the particles and set weights for full f
384 !!$ call self%sample( particle_group, params, xmin, Lx )
385 !!$ end if
386 !!$
387 !!$ counter = 0
388 !!$ ! Fill wi(2) with value of initial distribution at initial positions (g0)
389 !!$ ! and wi(3) with the delta f weights
390 !!$ do i_part = 1, particle_group%n_particles
391 !!$ q = particle_group%species%q
392 !!$ m = particle_group%species%m
393 !!$ xi = particle_group%get_x( i_part )
394 !!$ vi = particle_group%get_v( i_part )
395 !!$ wi = particle_group%get_weights( i_part )
396 !!$ ! TODO: Distinguish here between different initial sampling distributions
397 !!$ if( self%xiprofile ) then
398 !!$ if( present(map)) then
399 !!$ Rx = xi
400 !!$ Rx(1) = map%get_x1(xi) + m*vi(2)/q
401 !!$ Rx(1) = (Rx(1) - xmin(1))/Lx(1)
402 !!$ if( self%inverse ) then
403 !!$ wi(1) = wi(1) - params%eval_v_density( vi(1:params%dims(2)), Rx(1:params%dims(1)), m=particle_group%species%m )/params%eval_v_density( vi(1:params%dims(2)), xi(1:params%dims(1)), m=particle_group%species%m )*map%volume!control_variate%update_df_weight( Rx, vi, time0, wi(1), wi(2) )
404 !!$ else
405 !!$ wi(1) = wi(1) - params%eval_v_density( vi(1:params%dims(2)), Rx(1:params%dims(1)), m=particle_group%species%m )/params%eval_v_density( vi(1:params%dims(2)), xi(1:params%dims(1)), m=particle_group%species%m )*map%jacobian(xi)
406 !!$ end if
407 !!$ else
408 !!$ x = (xi - xmin)/Lx
409 !!$ Rx = x
410 !!$ Rx(1) = xi(1) + m*vi(2)/q
411 !!$ Rx(1) = (Rx(1) - xmin(1))/Lx(1)
412 !!$ wi(1) = wi(1) - params%eval_v_density( vi(1:params%dims(2)), Rx(1:params%dims(1)), m=particle_group%species%m )/params%eval_v_density( vi(1:params%dims(2)), x(1:params%dims(1)), m=particle_group%species%m )*product(Lx)!control_variate%update_df_weight( Rx, vi, time0, wi(1), wi(2) )
413 !!$ end if
414 !!$ else
415 !!$ if(delta_perturb) then
416 !!$ if( particle_group%species%q < 0._f64) then
417 !!$ if( x(1) > 0.5_f64-self%eps(1) .and. x(1) < 0.5_f64+self%eps(1) )then
418 !!$ if( x(2) > 0.5_f64-self%eps(2) .and. x(2) < 0.5_f64+self%eps(2) ) then
419 !!$ if(xi(3) > 0.5_f64-self%eps(3) .and. x(3) < 0.5_f64+self%eps(3) )then
420 !!$ vi = vi + self%a
421 !!$ counter = counter + 1
422 !!$ end if
423 !!$ end if
424 !!$ end if
425 !!$ end if
426 !!$ end if
427 !!$ if( present(map)) then
428 !!$ wi(1) = wi(1) - map%jacobian(xi)!control_variate%update_df_weight( xi, vi, time0, wi(1), wi(2) )
429 !!$ else
430 !!$ wi(1) = wi(1) - product(Lx)
431 !!$ end if
432 !!$ call particle_group%set_v(i_part, vi)
433 !!$ end if
434 !!$
435 !!$ call particle_group%set_weights( i_part, wi )
436 !!$ end do
437 !!$ if(delta_perturb) then
438 !!$ print*, 'delta_perturb_cv: ', self%eps, self%a, 'shifted particles:', counter
439 !!$ end if
440 !!$ end subroutine sample_lindf_particle_sampling
441 
443  subroutine sample_particle_sampling( self, particle_group, params, xmin, Lx, map )
444  class( sll_t_particle_sampling ), intent( inout ) :: self
445  class(sll_c_particle_group_base), target, intent(inout) :: particle_group
446  class( sll_c_distribution_params ), target, intent( inout ) :: params
447  sll_real64, intent(in) :: xmin(:)
448  sll_real64, intent(in) :: lx(:)
449  type(sll_t_mapping_3d), optional :: map
450 
451  !select type( params )
452  !type is( sll_t_params_cos_gaussian)
453 
454  if( present(map)) then
455 
456  if( self%symmetric == 0 ) then
457  call sample_particle_sampling_all_trafo( self, particle_group, params, xmin, lx, map )
458  elseif ( self%symmetric == 1 ) then
459  if ( params%dims(1) == 1 .and. params%dims(2) == 2 ) then
460  call sample_particle_sampling_sym_1d2v_trafo( self, particle_group, params, lx, map )
461  elseif ( params%dims(1) == 3 .and. params%dims(2) == 3 ) then
462  call sample_particle_sampling_sym_3d3v_trafo( self, particle_group, params, xmin, lx, map )
463  else
464  sll_error("sample_particle_sampling", "symmetric sampling not implemented for given dimension")
465  end if
466  else
467  sll_error("sample_particle_sampling", "this symmetric sampling is not implemented")
468  end if
469 
470  else
471  if( self%symmetric == 0 ) then
472  call sample_particle_sampling_all( self, particle_group, params, xmin, lx )
473  elseif ( self%symmetric == 1 ) then
474  if ( params%dims(1) == 1 .and. params%dims(2) == 2 ) then
475  call sample_particle_sampling_sym_1d2v( self, particle_group, params, xmin, lx )
476  elseif ( params%dims(1) == 3 .and. params%dims(2) == 3 ) then
477  if( self%uniform ) then
478  call sample_particle_sampling_sym_uni_3d3v( self, particle_group, params, xmin, lx )
479  else
480  call sample_particle_sampling_sym_3d3v( self, particle_group, params, xmin, lx )
481  end if
482  else
483  sll_error("sample_particle_sampling", "symmetric sampling not implemented for given dimension")
484  end if
485  elseif ( self%symmetric == 2) then
486  if ( params%dims(1) == 1 .and. params%dims(2) == 2 ) then
487  call sample_particle_sampling_sym2_1d2v( self, particle_group, params, xmin, lx )
488  else
489  sll_error("sample_particle_sampling", "symmetric sampling not implemented for given dimension")
490  end if
491  elseif( self%symmetric == 3) then
492  if ( params%dims(1) == 1 .and. params%dims(2) == 2 ) then
493  call sample_particle_sampling_unix_1d2v( self, particle_group, params, xmin, lx )
494  end if
495  end if
496  end if
497  !end select
498 
499  end subroutine sample_particle_sampling
500 
501 
503  subroutine sample_particle_sampling_all( self, particle_group, params, xmin, Lx )
504  type( sll_t_particle_sampling ), intent( inout ) :: self
505  class(sll_c_particle_group_base), target, intent( inout ) :: particle_group
506  class( sll_c_distribution_params ), target, intent( inout ) :: params
507  sll_real64, intent( in ) :: xmin(:)
508  sll_real64, intent( in ) :: lx(:)
509  !local variables
510  sll_int32 :: n_rnds
511  sll_real64 :: x(3), xi(3), v(3)
512  sll_int32 :: i_part
513  sll_int32 :: i_v, i_gauss
514  sll_real64, allocatable :: rdn(:)
515  sll_real64 :: wi(particle_group%n_weights)
516  sll_real64 :: rnd_no
517  sll_real64 :: delta(params%n_gaussians)
518  sll_real64 :: common_weight
519  sll_real64 :: r, rx(3)
520 
521  n_rnds = 0
522  if ( params%n_gaussians > 1 ) then
523  n_rnds = 1
524  end if
525 
526  do i_v=1,params%n_gaussians
527  delta(i_v) = sum(params%delta(1:i_v))
528  end do
529 
530 
531  n_rnds = n_rnds+params%dims(1)+params%dims(2)
532  allocate( rdn(1:params%dims(1)+params%dims(2)+1) )
533  rdn = 0.0_f64
534 
535  ! 1/Np in common weight
536  common_weight = particle_group%get_common_weight()
537  call particle_group%set_common_weight &
538  (common_weight/real(particle_group%n_total_particles, f64))
539 
540  if ( self%random_numbers == sll_p_random_numbers ) then
541  call random_seed(put=self%random_seed)
542  end if
543 
544  do i_part = 1, particle_group%n_particles
545  ! Generate Random or Sobol numbers on [0,1]
546  select case( self%random_numbers )
547  case( sll_p_sobol_numbers )
548  call sll_s_i8_sobol( int(n_rnds,8), self%sobol_seed, rdn(1:n_rnds))
549  case( sll_p_random_numbers )
550  call random_number( rdn(1:n_rnds) )
551  end select
552 
553  ! Transform rdn to the interval
554  xi(1:params%dims(1)) = rdn(1:params%dims(1))
555  x(1:params%dims(1)) = xmin + lx * xi(1:params%dims(1))
556 
557  ! Maxwellian distribution of the temperature
558  do i_v = 1,params%dims(2)
559  call sll_s_normal_cdf_inv( rdn(i_v+params%dims(1)), 0.0_f64, 1.0_f64, &
560  v(i_v))
561  end do
562  ! For multiple Gaussian, draw which one to take
563  rnd_no = rdn(params%dims(1)+params%dims(2)+1)
564  i_gauss = 1
565  do while( rnd_no > delta(i_gauss) )
566  i_gauss = i_gauss+1
567  end do
568  if( self%xiprofile ) then
569  select type(p=>params)
570  type is(sll_t_params_cos_gaussian_screwpinch)
571  if( particle_group%species%q > 0._f64) then
572  v(1:p%dims(2)) = v(1:p%dims(2)) * sqrt(p%profile%T_i(xi(1))/particle_group%species%m) + p%v_mean(:,i_gauss)
573  else
574  v(1:p%dims(2)) = v(1:p%dims(2)) * sqrt(p%profile%T_e(xi(1))/particle_group%species%m) + p%v_mean(:,i_gauss)
575  end if
576 
577  rx = xi
578  rx(1) = x(1) + particle_group%species%m*v(2)/particle_group%species%q
579  rx(1) = (rx(1) - xmin(1))/lx(1)
580 
581  !special perturbation in gyrocoordinates
582  wi(1) = p%eval_x_density(xi(1:p%dims(1)), v(1:p%dims(2)))*&
583  p%eval_v_density(v(1:p%dims(2)), rx(1:p%dims(1)), particle_group%species%m)/p%eval_v_density(v(1:p%dims(2)), xi(1:p%dims(1)), particle_group%species%m) *product(lx)
584  type is(sll_t_params_noise_gaussian)
585  if( particle_group%species%q > 0._f64) then
586  v(1:p%dims(2)) = v(1:p%dims(2)) * sqrt(p%profile%T_i(xi(1))/particle_group%species%m) + p%v_mean(:,i_gauss)
587  else
588  v(1:p%dims(2)) = v(1:p%dims(2)) * sqrt(p%profile%T_e(xi(1))/particle_group%species%m) + p%v_mean(:,i_gauss)
589  end if
590  wi(1) = p%eval_x_density(xi(1:p%dims(1))) *product(lx)
591  end select
592  else
593  v(1:params%dims(2)) = v(1:params%dims(2)) * params%v_thermal(:,i_gauss) + params%v_mean(:,i_gauss)
594  ! Set weight according to value of perturbation
595  wi(1) = params%eval_x_density(x(1:params%dims(1)),v(1:params%dims(2)))*product(lx)
596  end if
597  ! Copy the generated numbers to the particle
598  call particle_group%set_x(i_part, x)
599  call particle_group%set_v(i_part, v)
600  ! Set weights.
601  call particle_group%set_weights(i_part, &
602  wi)
603  end do
604 
605  if( particle_group%species%q < 0._f64) then
606  if( self%delta_perturb ) then
607  call delta_function_perturbation( self, particle_group, xmin, lx )
608  end if
609  end if
610 
611  end subroutine sample_particle_sampling_all
612 
613 
615  subroutine sample_particle_sampling_sym_1d2v( self, particle_group, params, xmin, Lx )
616  type( sll_t_particle_sampling ), intent( inout ) :: self
617  class(sll_c_particle_group_base), target, intent( inout ) :: particle_group
618  class( sll_c_distribution_params ), target, intent( in ) :: params
619  sll_real64, intent(in) :: xmin(:)
620  sll_real64, intent(in) :: lx(:)
621  !local variables
622  sll_int32 :: n_rnds
623  sll_real64 :: x(3),v(3)
624  sll_int32 :: i_part
625  sll_int32 :: i_v
626  sll_real64, allocatable :: rdn(:)
627  sll_real64 :: wi(1)
628  sll_real64 :: rnd_no
629  sll_int32 :: ip, i_gauss, counter
630  sll_real64 :: delta(params%n_gaussians)
631 
632  n_rnds = 0
633  if ( params%n_gaussians > 1 ) then
634  n_rnds = 1
635  end if
636 
637  do i_v=1,params%n_gaussians
638  delta(i_v) = sum(params%delta(1:i_v))
639  end do
640 
641  n_rnds = n_rnds+params%dims(1)+params%dims(2)
642  allocate( rdn(1:params%dims(1)+params%dims(2)+1) )
643  rdn = 0.0_f64
644 
645  ! 1/Np in common weight
646  call particle_group%set_common_weight &
647  (1.0_f64/real(particle_group%n_total_particles, f64))
648 
649  if ( self%random_numbers == sll_p_random_numbers ) then
650  call random_seed(put=self%random_seed)
651  end if
652 
653  do i_part = 1, particle_group%n_particles
654  ip = modulo(i_part, 8 )
655  if ( ip == 1) then
656  ! Generate Random or Sobol numbers on [0,1]
657  select case( self%random_numbers )
658  case( sll_p_sobol_numbers )
659  call sll_s_i8_sobol( int(n_rnds,8), self%sobol_seed, rdn(1:n_rnds))
660  case( sll_p_random_numbers )
661  call random_number( rdn(1:n_rnds) )
662  end select
663 
664  ! Transform rdn to the interval
665  x(1:params%dims(1)) = xmin + lx * rdn(1:params%dims(1))
666 
667  ! Set weight according to value of perturbation
668  wi(1) = params%eval_x_density(x(1:params%dims(1)))*product(lx)
669 
670  ! Maxwellian distribution of the temperature
671  do i_v = 1,params%dims(2)
672  call sll_s_normal_cdf_inv( rdn(i_v+params%dims(1)), 0.0_f64, 1.0_f64, &
673  v(i_v))
674  end do
675  ! For multiple Gaussian, draw which one to take
676  rnd_no = rdn(params%dims(1)+params%dims(2)+1)
677  i_gauss = 1
678  do while( rnd_no > delta(i_gauss) )
679  i_gauss = i_gauss+1
680  end do
681  v(1:params%dims(2)) = v(1:params%dims(2)) * params%v_thermal(:,i_gauss) + params%v_mean(:,i_gauss)
682  elseif ( ip == 5 ) then
683  x(1) = lx(1) - x(1) + 2.0_f64*xmin(1)
684  ! Set weight according to value of perturbation
685  wi(1) = params%eval_x_density(x(1:params%dims(1)))*product(lx)
686 
687  elseif ( modulo(ip,2) == 0 ) then
688  v(1) = -v(1) + 2.0_f64*params%v_mean(1,i_gauss)
689  else
690  v(2) = -v(2) + 2.0_f64*params%v_mean(2,i_gauss)
691  end if
692 
693  ! Copy the generated numbers to the particle
694  call particle_group%set_x(i_part, x)
695  call particle_group%set_v(i_part, v)
696  ! Set weights.
697  call particle_group%set_weights(i_part, &
698  wi)
699 
700  end do
701 
702  if( particle_group%species%q < 0._f64) then
703  if( self%delta_perturb ) then
704  call delta_function_perturbation( self, particle_group, xmin, lx )
705  end if
706  end if
707 
708  end subroutine sample_particle_sampling_sym_1d2v
709 
710 
711 
713  subroutine sample_particle_sampling_sym2_1d2v( self, particle_group, params, xmin, Lx )
714  type( sll_t_particle_sampling ), intent( inout ) :: self
715  class(sll_c_particle_group_base), target, intent(inout) :: particle_group
716  class( sll_c_distribution_params ), target, intent( inout ) :: params
717  sll_real64, intent(in) :: xmin(:)
718  sll_real64, intent(in) :: lx(:)
719  !local variables
720  sll_int32 :: n_rnds
721  sll_real64 :: x(3),v(3), xi(3)
722  sll_int32 :: i_part
723  sll_int32 :: i_v
724  sll_real64, allocatable :: rdn(:)
725  sll_real64 :: wi(1)
726  sll_real64 :: rnd_no
727  sll_int32 :: ip, i_gauss
728  sll_real64 :: delta(params%n_gaussians)
729  sll_real64 :: common_weight
730 
731  print*, 'new sampling'
732 
733  n_rnds = 0
734  if ( params%n_gaussians > 1 ) then
735  n_rnds = 1
736  end if
737 
738  do i_v=1,params%n_gaussians
739  delta(i_v) = sum(params%delta(1:i_v))
740  end do
741 
742  n_rnds = n_rnds+params%dims(1)+params%dims(2)
743  allocate( rdn(1:params%dims(1)+params%dims(2)+1) )
744  rdn = 0.0_f64
745 
746  ! 1/Np in common weight
747  common_weight = particle_group%get_common_weight()
748  call particle_group%set_common_weight &
749  (common_weight/real(particle_group%n_total_particles, f64))
750 
751  if ( self%random_numbers == sll_p_random_numbers ) then
752  call random_seed(put=self%random_seed)
753  end if
754 
755  do i_part = 1, particle_group%n_particles
756  ip = modulo(i_part, 2 )
757  if ( ip == 1) then
758  ! Generate Random or Sobol numbers on [0,1]
759  select case( self%random_numbers )
760  case( sll_p_sobol_numbers )
761  call sll_s_i8_sobol( int(n_rnds,8), self%sobol_seed, rdn(1:n_rnds))
762  case( sll_p_random_numbers )
763  call random_number( rdn(1:n_rnds) )
764  end select
765 
766  ! Transform rdn to the interval
767  xi(1:params%dims(1)) = rdn(1:params%dims(1))
768  x(1:params%dims(1)) = xmin + lx * xi(1:params%dims(1))
769 
770  ! Set weight according to value of perturbation
771  wi(1) = params%eval_x_density(x(1:params%dims(1)))*product(lx)
772 
773  ! Maxwellian distribution of the temperature
774  do i_v = 1,params%dims(2)
775  call sll_s_normal_cdf_inv( rdn(i_v+params%dims(1)), 0.0_f64, 1.0_f64, &
776  v(i_v))
777  end do
778  ! For multiple Gaussian, draw which one to take
779  rnd_no = rdn(params%dims(1)+params%dims(2)+1)
780  i_gauss = 1
781  do while( rnd_no > delta(i_gauss) )
782  i_gauss = i_gauss+1
783  end do
784  v(1:params%dims(2)) = v(1:params%dims(2)) * params%v_thermal(:,i_gauss) + params%v_mean(:,i_gauss)
785  else
786  x(1) = lx(1) - x(1) + 2.0_f64*xmin(1)
787  v(1) = -v(1) + 2.0_f64*params%v_mean(1,i_gauss)
788  v(2) = -v(2) + 2.0_f64*params%v_mean(2,i_gauss)
789  end if
790 
791  ! Copy the generated numbers to the particle
792  call particle_group%set_x(i_part, x)
793  call particle_group%set_v(i_part, v)
794  ! Set weights.
795  call particle_group%set_weights(i_part, &
796  wi)
797 
798  end do
799 
800  end subroutine sample_particle_sampling_sym2_1d2v
801 
803  subroutine sample_particle_sampling_unix( self, particle_group, params, xmin, Lx )
804  type( sll_t_particle_sampling ), intent( inout ) :: self
805  class(sll_c_particle_group_base), target, intent(inout) :: particle_group
806  class( sll_c_distribution_params ), target, intent( inout ) :: params
807  sll_real64, intent(in) :: xmin(:)
808  sll_real64, intent(in) :: lx(:)
809  !local variables
810  sll_int32 :: n_rnds
811  sll_real64 :: x(3),v(3), xi(3)
812  sll_int32 :: i_part
813  sll_int32 :: i_v
814  sll_real64, allocatable :: rdn(:)
815  sll_real64 :: wi(1)
816  sll_real64 :: rnd_no
817  sll_int32 :: i_gauss
818  sll_real64 :: delta(params%n_gaussians)
819  sll_real64 :: dx(params%dims(1))
820  sll_real64 :: common_weight
821 
822  print*, 'new sampling 2'
823 
824  n_rnds = 0
825  if ( params%n_gaussians > 1 ) then
826  n_rnds = 1
827  end if
828 
829  do i_v=1,params%n_gaussians
830  delta(i_v) = sum(params%delta(1:i_v))
831  end do
832 
833  n_rnds = n_rnds+params%dims(2)
834  allocate( rdn(1:params%dims(1)+params%dims(2)+1) )
835  rdn = 0.0_f64
836 
837  ! 1/Np in common weight
838  common_weight = particle_group%get_common_weight()
839  call particle_group%set_common_weight &
840  (common_weight/real(particle_group%n_total_particles, f64))
841 
842  if ( self%random_numbers == sll_p_random_numbers ) then
843  call random_seed(put=self%random_seed)
844  end if
845 
846  dx = lx/real(particle_group%n_particles/2,f64)
847 
848  do i_part = 1, particle_group%n_particles/4
849 
850  ! Set x value
851  x(1:params%dims(1)) = xmin + real(i_part-1, f64)*dx
852  ! Set weight according to value of perturbation
853  wi(1) = params%eval_x_density(x(1:params%dims(1)))*product(lx)
854 
855  ! Generate first draw for v
856  ! Generate Random or Sobol numbers on [0,1]
857  select case( self%random_numbers )
858  case( sll_p_sobol_numbers )
859  call sll_s_i8_sobol( int(n_rnds,8), self%sobol_seed, rdn(1:n_rnds))
860  case( sll_p_random_numbers )
861  call random_number( rdn(1:n_rnds) )
862  end select
863  ! Maxwellian distribution of the temperature
864  do i_v = 1,params%dims(2)
865  call sll_s_normal_cdf_inv( rdn(i_v), 0.0_f64, 1.0_f64, &
866  v(i_v))
867  end do
868  ! For multiple Gaussian, draw which one to take
869  rnd_no = rdn(params%dims(2)+1)
870  i_gauss = 1
871  do while( rnd_no > delta(i_gauss) )
872  i_gauss = i_gauss+1
873  end do
874  v(1:params%dims(2)) = v(1:params%dims(2)) * params%v_thermal(:,i_gauss) + params%v_mean(:,i_gauss)
875 
876  ! Copy the generated numbers to the particle
877  call particle_group%set_x(2*i_part-1, x)
878  call particle_group%set_v(2*i_part-1, v)
879  ! Set weights.
880  call particle_group%set_weights(2*i_part-1, &
881  wi)
882 
883  ! Generate second random number
884  ! Generate Random or Sobol numbers on [0,1]
885  select case( self%random_numbers )
886  case( sll_p_sobol_numbers )
887  call sll_s_i8_sobol( int(n_rnds,8), self%sobol_seed, rdn(1:n_rnds))
888  case( sll_p_random_numbers )
889  call random_number( rdn(1:n_rnds) )
890  end select
891  ! Maxwellian distribution of the temperature
892  do i_v = 1,params%dims(2)
893  call sll_s_normal_cdf_inv( rdn(i_v), 0.0_f64, 1.0_f64, &
894  v(i_v))
895  end do
896  ! For multiple Gaussian, draw which one to take
897  rnd_no = rdn(params%dims(2)+1)
898  i_gauss = 1
899  do while( rnd_no > delta(i_gauss) )
900  i_gauss = i_gauss+1
901  end do
902  v(1:params%dims(2)) = v(1:params%dims(2)) * params%v_thermal(:,i_gauss) + params%v_mean(:,i_gauss)
903 
904  ! Copy the generated numbers to the particle
905  call particle_group%set_x(2*i_part, x)
906  call particle_group%set_v(2*i_part, v)
907  ! Set weights.
908  call particle_group%set_weights(2*i_part, &
909  wi)
910  end do
911 
912  x = xmin + lx*0.5_f64
913  call particle_group%set_x(1, x)
914 
915  ! Add the negative particles
916  do i_part = 1,particle_group%n_particles/2
917 
918  x = particle_group%get_x( i_part )
919  v = -particle_group%get_v( i_part )
920 
921  x(1:params%dims(1)) = lx - x(1:params%dims(1)) + 2.0_f64*xmin
922  ! xi(1:params%dims(1)) = 1._f64 - xi(1:params%dims(1))
923  wi(1) = params%eval_x_density(x(1:params%dims(1)))*product(lx)
924 
925  ! Copy the generated numbers to the particle
926  call particle_group%set_x(i_part+ particle_group%n_particles/2, x)
927  call particle_group%set_v(i_part+ particle_group%n_particles/2, v)
928  ! Set weights.
929  call particle_group%set_weights(i_part+ particle_group%n_particles/2, &
930  wi)
931 
932  end do
933 
934  end subroutine sample_particle_sampling_unix
935 
936 
938  subroutine sample_particle_sampling_unix_1d2v( self, particle_group, params, xmin, Lx )
939  type( sll_t_particle_sampling ), intent( inout ) :: self
940  class(sll_c_particle_group_base), target, intent(inout) :: particle_group
941  class( sll_c_distribution_params ), target, intent( inout ) :: params
942  sll_real64, intent(in) :: xmin(:)
943  sll_real64, intent(in) :: lx(:)
944  !local variables
945  sll_int32 :: n_rnds
946  sll_real64 :: x(3),v(3), xi(3)
947  sll_int32 :: i_part
948  sll_int32 :: j
949  sll_int32 :: i_v
950  sll_real64, allocatable :: rdn(:)
951  sll_real64 :: wi(1)
952  sll_real64 :: rnd_no
953  sll_int32 :: i_gauss
954  sll_real64 :: delta(params%n_gaussians)
955  sll_real64 :: dx(params%dims(1))
956  sll_real64 :: common_weight
957 
958  print*, 'new sampling 2'
959 
960  n_rnds = 0
961  if ( params%n_gaussians > 1 ) then
962  n_rnds = 1
963  end if
964 
965  do i_v=1,params%n_gaussians
966  delta(i_v) = sum(params%delta(1:i_v))
967  end do
968 
969  n_rnds = n_rnds+params%dims(2)
970  allocate( rdn(1:params%dims(1)+params%dims(2)+1) )
971  rdn = 0.0_f64
972 
973  ! 1/Np in common weight
974  common_weight = particle_group%get_common_weight()
975  call particle_group%set_common_weight &
976  (common_weight/real(particle_group%n_total_particles, f64))
977 
978  if ( self%random_numbers == sll_p_random_numbers ) then
979  call random_seed(put=self%random_seed)
980  end if
981 
982  dx = lx/real(particle_group%n_particles/4,f64)
983 
984  do i_part = 1, particle_group%n_particles/8+1
985 
986  ! Set x value
987  x(1:params%dims(1)) = xmin + real(i_part-1,f64)*dx
988  ! Set weight according to value of perturbation
989  wi(1) = params%eval_x_density(x(1:params%dims(1)))*product(lx)
990 
991  ! Generate first draw for v
992  ! Generate Random or Sobol numbers on [0,1]
993  select case( self%random_numbers )
994  case( sll_p_sobol_numbers )
995  call sll_s_i8_sobol( int(n_rnds,8), self%sobol_seed, rdn(1:n_rnds))
996  case( sll_p_random_numbers )
997  call random_number( rdn(1:n_rnds) )
998  end select
999  ! Maxwellian distribution of the temperature
1000  do i_v = 1,params%dims(2)
1001  call sll_s_normal_cdf_inv( rdn(i_v), 0.0_f64, 1.0_f64, &
1002  v(i_v))
1003  end do
1004  ! For multiple Gaussian, draw which one to take
1005  rnd_no = rdn(params%dims(2)+1)
1006  i_gauss = 1
1007  do while( rnd_no > delta(i_gauss) )
1008  i_gauss = i_gauss+1
1009  end do
1010  v(1:params%dims(2)) = v(1:params%dims(2)) * params%v_thermal(:,i_gauss) + params%v_mean(:,i_gauss)
1011 
1012  ! Copy the generated numbers to the particle
1013  do j = 0,3
1014  call particle_group%set_x(4*i_part-j, x)
1015  ! Set weights.
1016  call particle_group%set_weights(4*i_part-j, &
1017  wi)
1018  end do
1019 
1020  call particle_group%set_v(4*i_part-3, v)
1021  call particle_group%set_v(4*i_part-2, [-v(1),v(2),v(3)])
1022  call particle_group%set_v(4*i_part-1, [v(1),-v(2),v(3)])
1023  call particle_group%set_v(4*i_part, [-v(1),-v(2),v(3)])
1024 
1025 
1026  end do
1027 
1028  ! Add the negative particles
1029  do i_part = 5,particle_group%n_particles/2
1030 
1031  x = particle_group%get_x( i_part )
1032  v = -particle_group%get_v( i_part )
1033 
1034  x(1:params%dims(1)) = lx - x(1:params%dims(1)) + 2.0_f64*xmin
1035  !xi(1:params%dims(1)) = 1._f64 - xi(1:params%dims(1))
1036  wi(1) = params%eval_x_density(x(1:params%dims(1)))*product(lx)
1037 
1038  ! Copy the generated numbers to the particle
1039  call particle_group%set_x(i_part+ particle_group%n_particles/2, x)
1040  call particle_group%set_v(i_part+ particle_group%n_particles/2, v)
1041  ! Set weights.
1042  call particle_group%set_weights(i_part+ particle_group%n_particles/2, &
1043  wi)
1044 
1045  end do
1046 
1047  end subroutine sample_particle_sampling_unix_1d2v
1048 
1049 
1050 
1052  subroutine sample_particle_sampling_sym_3d3v( self, particle_group, params, xmin, Lx )
1053  type( sll_t_particle_sampling ), intent( inout ) :: self
1054  class(sll_c_particle_group_base), target, intent(inout) :: particle_group
1055  class( sll_c_distribution_params ), target, intent( in ) :: params
1056  sll_real64, intent(in) :: xmin(:)
1057  sll_real64, intent(in) :: lx(:)
1058  !local variables
1059  sll_int32 :: n_rnds
1060  sll_real64 :: x(3),v(3)
1061  sll_int32 :: i_part
1062  sll_int32 :: i_v
1063  sll_real64, allocatable :: rdn(:)
1064  sll_real64 :: wi(1)
1065  sll_real64 :: rnd_no
1066  sll_int32 :: ip, i_gauss
1067  sll_real64 :: delta(params%n_gaussians)
1068 
1069  n_rnds = 0
1070  if ( params%n_gaussians > 1 ) then
1071  n_rnds = 1
1072  end if
1073 
1074  do i_v=1,params%n_gaussians
1075  delta(i_v) = sum(params%delta(1:i_v))
1076  end do
1077 
1078  n_rnds = n_rnds+params%dims(1)+params%dims(2)
1079  allocate( rdn(1:params%dims(1)+params%dims(2)+1) )
1080  rdn = 0.0_f64
1081 
1082  ! 1/Np in common weight
1083  call particle_group%set_common_weight &
1084  (1.0_f64/real(particle_group%n_total_particles, f64))
1085 
1086  if ( self%random_numbers == sll_p_random_numbers ) then
1087  call random_seed(put=self%random_seed)
1088  end if
1089 
1090  do i_part = 1, particle_group%n_particles
1091  ip = modulo(i_part, 64 )
1092  if ( ip == 1) then
1093  ! Generate Random or Sobol numbers on [0,1]
1094  select case( self%random_numbers )
1095  case( sll_p_sobol_numbers )
1096  call sll_s_i8_sobol( int(n_rnds,8), self%sobol_seed, rdn(1:n_rnds))
1097  case( sll_p_random_numbers )
1098  call random_number( rdn(1:n_rnds) )
1099  end select
1100 
1101  ! Transform rdn to the interval
1102  x(1:params%dims(1)) = xmin + lx * rdn(1:params%dims(1))
1103 
1104  ! Set weight according to value of perturbation
1105  wi(1) = params%eval_x_density(x(1:params%dims(1)))*product(lx)
1106 
1107  ! Maxwellian distribution of the temperature
1108  do i_v = 1,params%dims(2)
1109  call sll_s_normal_cdf_inv( rdn(i_v+params%dims(1)), 0.0_f64, 1.0_f64, &
1110  v(i_v))
1111  end do
1112  ! For multiple Gaussian, draw which one to take
1113  rnd_no = rdn(params%dims(1)+params%dims(2)+1)
1114  i_gauss = 1
1115  do while( rnd_no > delta(i_gauss) )
1116  i_gauss = i_gauss+1
1117  end do
1118  v(1:params%dims(2)) = v(1:params%dims(2)) * params%v_thermal(:,i_gauss) + params%v_mean(:,i_gauss)
1119  elseif ( modulo(i_part, 2) == 0 ) then
1120  v(3) = - v(3) + 2.0_f64*params%v_mean(3,i_gauss)
1121  elseif ( modulo(i_part, 4) == 3 ) then
1122  v(2) = - v(2) + 2.0_f64*params%v_mean(2,i_gauss)
1123  elseif ( modulo(i_part, 8) == 5 ) then
1124  v(1) = - v(1) + 2.0_f64*params%v_mean(1,i_gauss)
1125  elseif ( modulo(i_part, 16) == 9 ) then
1126  x(3) = lx(3) - x(3) + 2.0_f64*xmin(3)
1127  ! Set weight according to value of perturbation
1128  wi(1) = params%eval_x_density(x(1:params%dims(1)))*product(lx)
1129  elseif ( modulo(i_part, 32) == 17 ) then
1130  x(2) = lx(2) - x(2) + 2.0_f64*xmin(2)
1131  ! Set weight according to value of perturbation
1132  wi(1) = params%eval_x_density(x(1:params%dims(1)))*product(lx)
1133  elseif ( modulo(i_part, 64) == 33 ) then
1134  x(1) = lx(1) - x(1) + 2.0_f64*xmin(1)
1135  ! Set weight according to value of perturbation
1136  wi(1) = params%eval_x_density(x(1:params%dims(1)))*product(lx)
1137  end if
1138 
1139  ! Copy the generated numbers to the particle
1140  call particle_group%set_x(i_part, x)
1141  call particle_group%set_v(i_part, v)
1142  ! Set weights.
1143  call particle_group%set_weights(i_part, &
1144  wi)
1145 
1146  end do
1147 
1148  if( particle_group%species%q < 0._f64) then
1149  if( self%delta_perturb ) then
1150  call delta_function_perturbation( self, particle_group, xmin, lx )
1151  end if
1152  end if
1153 
1154  end subroutine sample_particle_sampling_sym_3d3v
1155 
1156 
1158  subroutine sample_particle_sampling_sym_uni_3d3v( self, particle_group, params, xmin, Lx )
1159  type( sll_t_particle_sampling ), intent( inout ) :: self
1160  class(sll_c_particle_group_base), target, intent(inout) :: particle_group
1161  class( sll_c_distribution_params ), target, intent( inout ) :: params
1162  sll_real64, intent(in) :: xmin(:)
1163  sll_real64, intent(in) :: lx(:)
1164  !local variables
1165  sll_int32 :: n_rnds
1166  sll_real64 :: x(3), v(3), xi(3)
1167  sll_int32 :: i_part, i_v
1168  sll_real64, allocatable :: rdn(:)
1169  sll_real64 :: wi(1)
1170  sll_real64 :: rnd_no
1171  sll_int32 :: ip, i_gauss
1172  sll_real64 :: v_min(3), lv(3)
1173 
1174  n_rnds = 0
1175  if ( params%n_gaussians > 1 ) then
1176  n_rnds = 1
1177  end if
1178 
1179 
1180  n_rnds = n_rnds+params%dims(1)+params%dims(2)
1181  allocate( rdn(1:params%dims(1)+params%dims(2)+1) )
1182  rdn = 0.0_f64
1183 
1184  do i_v = 1, 3
1185  v_min(i_v) = minval(params%v_mean(i_v,:)) -3._f64 * maxval(params%v_thermal(i_v,:))
1186  lv(i_v) = v_min(i_v) - maxval(params%v_mean(i_v,:)) +3._f64 * maxval(params%v_thermal(i_v,:))
1187  end do
1188 
1189  ! 1/Np in common weight
1190  call particle_group%set_common_weight &
1191  (1.0_f64/real(particle_group%n_total_particles, f64))
1192 
1193  if ( self%random_numbers == sll_p_random_numbers ) then
1194  call random_seed(put=self%random_seed)
1195  end if
1196 
1197  do i_part = 1, particle_group%n_particles
1198  ip = modulo(i_part, 64 )
1199  if ( ip == 1) then
1200  ! Generate Random or Sobol numbers on [0,1]
1201  select case( self%random_numbers )
1202  case( sll_p_sobol_numbers )
1203  call sll_s_i8_sobol( int(n_rnds,8), self%sobol_seed, rdn(1:n_rnds))
1204  case( sll_p_random_numbers )
1205  call random_number( rdn(1:n_rnds) )
1206  end select
1207 
1208  ! Transform rdn to the interval, sample x in [x_min, x_min+L_x]
1209  xi(1:params%dims(1)) = rdn(1:params%dims(1))
1210  x(1:params%dims(1)) = xmin + lx * xi(1:params%dims(1))
1211 
1212  ! Sample v uniform in [v_min, v_min+Lv]
1213  v(1:params%dims(2)) = v_min + lv * rdn(params%dims(1)+1:params%dims(1)+params%dims(2))
1214 
1215  ! Set weight according to value of perturbation
1216  wi(1) = params%eval_xv_density(x(1:params%dims(1)), v(1:params%dims(2)))*product(lx)*product(lv)
1217  elseif ( modulo(i_part, 2) == 0 ) then
1218  v(3) = - v(3) + 2.0_f64*params%v_mean(3,i_gauss)
1219  elseif ( modulo(i_part, 4) == 3 ) then
1220  v(2) = - v(2) + 2.0_f64*params%v_mean(2,i_gauss)
1221  elseif ( modulo(i_part, 8) == 5 ) then
1222  v(1) = - v(1) + 2.0_f64*params%v_mean(1,i_gauss)
1223  elseif ( modulo(i_part, 16) == 9 ) then
1224  x(3) = lx(3) - x(3) + 2.0_f64*xmin(3)
1225 
1226  ! Set weight according to value of perturbation
1227  wi(1) = params%eval_x_density(x(1:params%dims(1)))*product(lx)
1228  elseif ( modulo(i_part, 32) == 17 ) then
1229  x(2) = lx(2) - x(2) + 2.0_f64*xmin(2)
1230 
1231  ! Set weight according to value of perturbation
1232  wi(1) = params%eval_x_density(x(1:params%dims(1)))*product(lx)
1233  elseif ( modulo(i_part, 64) == 33 ) then
1234  x(1) = lx(1) - x(1) + 2.0_f64*xmin(1)
1235 
1236  ! Set weight according to value of perturbation
1237  wi(1) = params%eval_x_density(x(1:params%dims(1)))*product(lx)
1238  end if
1239 
1240  ! Copy the generated numbers to the particle
1241  call particle_group%set_x(i_part, x)
1242  call particle_group%set_v(i_part, v)
1243  ! Set weights.
1244  call particle_group%set_weights(i_part, &
1245  wi)
1246 
1247  end do
1248 
1249  end subroutine sample_particle_sampling_sym_uni_3d3v
1250 
1251 
1253  subroutine sample_particle_sampling_all_trafo( self, particle_group, params, xmin, Lx, map )
1254  type(sll_t_particle_sampling), intent( inout ) :: self
1255  class(sll_c_particle_group_base), target, intent( inout ) :: particle_group
1256  class(sll_c_distribution_params), target, intent( inout ) :: params
1257  sll_real64, intent( in ) :: xmin(:)
1258  sll_real64, intent( in ) :: lx(:)
1259  type(sll_t_mapping_3d), intent( inout ) :: map
1260  !local variables
1261  sll_int32 :: n_rnds
1262  sll_real64 :: x(3), xi(3), v(3)
1263  sll_int32 :: i_part
1264  sll_int32 :: i_v, i_gauss
1265  sll_real64, allocatable :: rdn(:)
1266  sll_real64 :: wi(particle_group%n_weights)
1267  sll_real64 :: rnd_no
1268  sll_real64 :: delta(params%n_gaussians)
1269  sll_real64 :: w_total(1), w_local(1)
1270  sll_real64 :: common_weight
1271  sll_real64 :: r, rx(3)
1272 
1273  w_local=0._f64
1274  w_total=0._f64
1275  n_rnds = 0
1276  if ( params%n_gaussians > 1 ) then
1277  n_rnds = 1
1278  end if
1279 
1280  do i_v=1,params%n_gaussians
1281  delta(i_v) = sum(params%delta(1:i_v))
1282  end do
1283 
1284  n_rnds = n_rnds+params%dims(1)+params%dims(2)
1285  allocate( rdn(1:params%dims(1)+params%dims(2)+1) )
1286  rdn = 0.0_f64
1287 
1288  ! 1/Np in common weight
1289  common_weight = particle_group%get_common_weight()
1290  call particle_group%set_common_weight &
1291  (common_weight/real(particle_group%n_total_particles, f64))
1292 
1293  if ( self%random_numbers == sll_p_random_numbers ) then
1294  call random_seed(put=self%random_seed)
1295  end if
1296 
1297  do i_part = 1, particle_group%n_particles
1298  ! Generate Random or Sobol numbers on [0,1]
1299  select case( self%random_numbers )
1300  case( sll_p_sobol_numbers )
1301  call sll_s_i8_sobol( int(n_rnds,8), self%sobol_seed, rdn(1:n_rnds))
1302  case( sll_p_random_numbers )
1303  call random_number( rdn(1:n_rnds) )
1304  end select
1305 
1306  xi = 0._f64
1307  x = 0._f64
1308  v = 0._f64
1309 
1310  if( self%inverse ) then
1311  ! Transform rdn to the interval
1312  x(1:2) = 2._f64*map%params(2) * rdn(1:2) - map%params(2)
1313  x(3) = map%params(3) * rdn(3)
1314  r = sqrt(x(1)**2 +x(2)**2)
1315  do while( r > map%params(2) .or. r < map%params(1) )
1316  select case( self%random_numbers )
1317  case( sll_p_sobol_numbers )
1318  call sll_s_i8_sobol( int(n_rnds,8), self%sobol_seed, rdn(1:n_rnds))
1319  case( sll_p_random_numbers )
1320  call random_number( rdn(1:n_rnds) )
1321  end select
1322  ! Transform rdn to the interval
1323  x(1:2) = 2._f64*map%params(2) * rdn(1:2) - map%params(2)
1324  x(3) = map%params(3) * rdn(3)
1325  ! inverse transformation of the uniformly sampled particles in physical coordinates
1326  r = sqrt(x(1)**2 +x(2)**2)
1327  end do
1328  xi = map%get_xi(x)
1329  ! inverse transformation of the uniformly sampled particles in physical coordinates
1330  ! Set weight according to value of perturbation
1331  wi(1) = params%eval_x_density(x) * map%volume
1332  else
1333  xi(1:params%dims(1)) = rdn(1:params%dims(1))
1334  x = map%get_x(xi)
1335  end if
1336  ! Maxwellian distribution of the temperature
1337  do i_v = 1,params%dims(2)
1338  call sll_s_normal_cdf_inv( rdn(i_v+params%dims(1)), 0.0_f64, 1.0_f64, &
1339  v(i_v))
1340  end do
1341  ! For multiple Gaussian, draw which one to take
1342  rnd_no = rdn(params%dims(1)+params%dims(2)+1)
1343  i_gauss = 1
1344  do while( rnd_no > delta(i_gauss) )
1345  i_gauss = i_gauss+1
1346  end do
1347  if( self%xiprofile ) then
1348  select type(p=>params)
1349  type is(sll_t_params_cos_gaussian_screwpinch)
1350  if( particle_group%species%q > 0._f64) then
1351  v(1:p%dims(2)) = v(1:p%dims(2)) * sqrt(p%profile%T_i(xi(1))/particle_group%species%m) + p%v_mean(:,i_gauss)
1352  else
1353  v(1:p%dims(2)) = v(1:p%dims(2)) * sqrt(p%profile%T_e(xi(1))/particle_group%species%m) + p%v_mean(:,i_gauss)
1354  end if
1355  rx = xi
1356  rx(1) = x(1) + particle_group%species%m*v(2)/particle_group%species%q
1357  rx(1) = (rx(1) - xmin(1))/lx(1)
1358 
1359  !special perturbation in gyrocoordinates
1360  wi(1) = p%eval_x_density(xi(1:p%dims(1)), v(1:p%dims(2)))*&
1361  p%eval_v_density(v(1:p%dims(2)), rx(1:p%dims(1)), particle_group%species%m)/p%eval_v_density(v(1:p%dims(2)), xi(1:p%dims(1)), particle_group%species%m) *map%jacobian(xi)
1362  class default
1363  wi(1) = params%eval_x_density( xi(1:params%dims(1)) )*map%jacobian(xi)
1364  end select
1365  else
1366  v(1:params%dims(2)) = v(1:params%dims(2)) * params%v_thermal(:,i_gauss) + params%v_mean(:,i_gauss)
1367  ! Set weight according to value of perturbation
1368  wi(1) = params%eval_x_density( x(1:params%dims(1)) )*map%jacobian(xi)
1369  end if
1370  ! Copy the generated numbers to the particle
1371  call particle_group%set_x(i_part, xi)
1372  call particle_group%set_v(i_part, v)
1373  ! Set weights.
1374  call particle_group%set_weights(i_part, &
1375  wi)
1376  !w_local=w_local+wi(1)
1377  end do
1378 
1379  if( particle_group%species%q < 0._f64) then
1380  if( self%delta_perturb ) then
1381  call delta_function_perturbation( self, particle_group, [0._f64,0._f64,0._f64], [1._f64,1._f64,1._f64] )
1382  end if
1383  end if
1384 
1385  end subroutine sample_particle_sampling_all_trafo
1386 
1387  subroutine delta_function_perturbation( self, particle_group, xmin, Lx )
1388  type(sll_t_particle_sampling), intent( inout ) :: self
1389  class(sll_c_particle_group_base), target, intent( inout ) :: particle_group
1390  sll_real64, intent(in) :: xmin(:)
1391  sll_real64, intent(in) :: lx(:)
1392  !local variables
1393  sll_int32 :: n_rnds
1394  sll_real64 :: xi(3), vi(3)
1395  sll_int32 :: i_part, counter
1396 
1397  counter = 0
1398  do i_part = 1, particle_group%n_particles
1399  xi = particle_group%get_x(i_part)
1400  vi = particle_group%get_v(i_part)
1401 
1402  xi = (xi - xmin)/lx
1403  if( xi(1) > 0.5_f64-self%eps(1) .and. xi(1) < 0.5_f64+self%eps(1) )then
1404  if( xi(2) > 0.5_f64-self%eps(2) .and. xi(2) < 0.5_f64+self%eps(2) ) then
1405  if(xi(3) > 0.5_f64-self%eps(3) .and. xi(3) < 0.5_f64+self%eps(3) )then
1406  vi = vi + self%a
1407  counter = counter + 1
1408  end if
1409  end if
1410  end if
1411  call particle_group%set_v(i_part, vi)
1412  end do
1413  print*, 'delta_perturb, ', self%eps, self%a, 'shifted particles:', counter
1414  end subroutine delta_function_perturbation
1415 
1416 
1418  subroutine sample_particle_sampling_sym_1d2v_trafo( self, particle_group, params, Lx, map )
1419  type(sll_t_particle_sampling), intent( inout ) :: self
1420  class(sll_c_particle_group_base), target, intent( inout ) :: particle_group
1421  class(sll_c_distribution_params), target, intent( inout ) :: params
1422  sll_real64, intent( in ) :: lx(:)
1423  type(sll_t_mapping_3d), intent( inout ) :: map
1424  !local variables
1425  sll_int32 :: n_rnds
1426  sll_real64 :: x(3),v(3), xi(3)
1427  sll_int32 :: i_part
1428  sll_int32 :: i_v
1429  sll_real64, allocatable :: rdn(:)
1430  sll_real64 :: wi(1)
1431  sll_real64 :: rnd_no
1432  sll_int32 :: ip, i_gauss
1433  sll_real64 :: delta(params%n_gaussians)
1434  sll_real64 :: common_weight
1435 
1436  n_rnds = 0
1437  if ( params%n_gaussians > 1 ) then
1438  n_rnds = 1
1439  end if
1440 
1441  do i_v=1,params%n_gaussians
1442  delta(i_v) = sum(params%delta(1:i_v))
1443  end do
1444 
1445  n_rnds = n_rnds+params%dims(1)+params%dims(2)
1446  allocate( rdn(1:params%dims(1)+params%dims(2)+1) )
1447  rdn = 0.0_f64
1448 
1449  ! 1/Np in common weight
1450  common_weight = particle_group%get_common_weight()
1451  call particle_group%set_common_weight &
1452  (common_weight/real(particle_group%n_total_particles, f64))
1453 
1454  if ( self%random_numbers == sll_p_random_numbers ) then
1455  call random_seed(put=self%random_seed)
1456  end if
1457 
1458  do i_part = 1, particle_group%n_particles
1459  ip = modulo(i_part, 8 )
1460  if ( ip == 1) then
1461  ! Generate Random or Sobol numbers on [0,1]
1462  select case( self%random_numbers )
1463  case( sll_p_sobol_numbers )
1464  call sll_s_i8_sobol( int(n_rnds,8), self%sobol_seed, rdn(1:n_rnds))
1465  case( sll_p_random_numbers )
1466  call random_number( rdn(1:n_rnds) )
1467  end select
1468 
1469  xi = 0._f64
1470  x = 0._f64
1471  v = 0._f64
1472 
1473  ! Transform rdn to the interval
1474  xi(1:params%dims(1)) = rdn(1:params%dims(1))
1475  x(1) = map%get_x1(xi)
1476  ! Set weight according to value of perturbation
1477  wi(1) = params%eval_x_density(x(1:params%dims(1)))*abs(map%jacobian(xi))
1478 
1479  ! Maxwellian distribution of the temperature
1480  do i_v = 1,params%dims(2)
1481  call sll_s_normal_cdf_inv( rdn(i_v+params%dims(1)), 0.0_f64, 1.0_f64, &
1482  v(i_v))
1483  end do
1484  ! For multiple Gaussian, draw which one to take
1485  rnd_no = rdn(params%dims(1)+params%dims(2)+1)
1486  i_gauss = 1
1487  do while( rnd_no > delta(i_gauss) )
1488  i_gauss = i_gauss+1
1489  end do
1490  v(1:params%dims(2)) = v(1:params%dims(2)) * params%v_thermal(:,i_gauss) + params%v_mean(:,i_gauss)
1491  elseif ( ip == 5 ) then
1492  xi(1) = 1._f64 - xi(1)
1493  x = 0._f64
1494  x(1) = map%get_x1(xi)
1495  ! Set weight according to value of perturbation
1496  wi(1) = params%eval_x_density(x(1:params%dims(1)))*abs(map%jacobian(xi))
1497  elseif ( modulo(ip,2) == 0 ) then
1498  v(1) = -v(1) + 2.0_f64*params%v_mean(1,i_gauss)
1499  else
1500  v(2) = -v(2) + 2.0_f64*params%v_mean(2,i_gauss)
1501  end if
1502 
1503  ! Copy the generated numbers to the particle
1504  call particle_group%set_x(i_part, xi)
1505  call particle_group%set_v(i_part, v)
1506  ! Set weights.
1507  call particle_group%set_weights(i_part, &
1508  wi)
1509  end do
1510 
1511  end subroutine sample_particle_sampling_sym_1d2v_trafo
1512 
1513 
1515  subroutine sample_particle_sampling_sym_3d3v_trafo( self, particle_group, params, xmin, Lx, map )
1516  type(sll_t_particle_sampling), intent( inout ) :: self
1517  class(sll_c_particle_group_base), target, intent( inout ) :: particle_group
1518  class(sll_c_distribution_params), target, intent( inout ) :: params
1519  sll_real64, intent( in ) :: xmin(:)
1520  sll_real64, intent( in ) :: lx(:)
1521  type(sll_t_mapping_3d), intent( inout ) :: map
1522  !local variables
1523  sll_int32 :: n_rnds
1524  sll_real64 :: x(3), xi(3), v(3)
1525  sll_int32 :: i_part
1526  sll_int32 :: i_v
1527  sll_real64, allocatable :: rdn(:)
1528  sll_real64 :: wi(particle_group%n_weights)
1529  sll_real64 :: rnd_no
1530  sll_int32 :: ip, i_gauss
1531  sll_real64 :: delta(params%n_gaussians)
1532  sll_real64 :: w_total(1), w_local(1)
1533  sll_real64 :: common_weight, weight
1534  sll_real64 :: r
1535 
1536  w_local=0._f64
1537  w_total=0._f64
1538  n_rnds = 0
1539  if ( params%n_gaussians > 1 ) then
1540  n_rnds = 1
1541  end if
1542 
1543  do i_v=1,params%n_gaussians
1544  delta(i_v) = sum(params%delta(1:i_v))
1545  end do
1546 
1547  n_rnds = n_rnds+params%dims(1)+params%dims(2)
1548  allocate( rdn(1:params%dims(1)+params%dims(2)+1) )
1549  rdn = 0.0_f64
1550 
1551  ! 1/Np in common weight
1552  common_weight = particle_group%get_common_weight()
1553  call particle_group%set_common_weight &
1554  (common_weight/real(particle_group%n_total_particles, f64))
1555 
1556  if ( self%random_numbers == sll_p_random_numbers ) then
1557  call random_seed(put=self%random_seed)
1558  end if
1559 
1560  do i_part = 1, particle_group%n_particles
1561  ip = modulo(i_part, 4 )
1562  if ( ip == 1) then
1563  ! Generate Random or Sobol numbers on [0,1]
1564  select case( self%random_numbers )
1565  case( sll_p_sobol_numbers )
1566  call sll_s_i8_sobol( int(n_rnds,8), self%sobol_seed, rdn(1:n_rnds))
1567  case( sll_p_random_numbers )
1568  call random_number( rdn(1:n_rnds) )
1569  end select
1570 
1571  if( self%inverse ) then
1572  ! Transform rdn to the interval
1573  x(1:params%dims(1)-1) = 2._f64*map%params(2) * rdn(1:params%dims(1)-1) - map%params(2)
1574  x(params%dims(1)) = map%params(3) * rdn(params%dims(1))
1575  ! inverse transformation of the uniformly sampled particles in physical coordinates
1576  r = sqrt(x(1)**2 +x(2)**2)
1577  do while( r > map%params(2) .or. r < map%params(1) )
1578  select case( self%random_numbers )
1579  case( sll_p_sobol_numbers )
1580  call sll_s_i8_sobol( int(n_rnds,8), self%sobol_seed, rdn(1:n_rnds))
1581  case( sll_p_random_numbers )
1582  call random_number( rdn(1:n_rnds) )
1583  end select
1584  ! Transform rdn to the interval
1585  x(1:params%dims(1)-1) = 2._f64*map%params(2) * rdn(1:params%dims(1)-1) - map%params(2)
1586  x(params%dims(1)) = map%params(3) * rdn(params%dims(1))
1587  ! inverse transformation of the uniformly sampled particles in physical coordinates
1588  r = sqrt(x(1)**2 +x(2)**2)
1589  end do
1590  xi = map%get_xi(x)
1591  ! Set weight according to value of perturbation
1592  wi(1) = params%eval_x_density(x) * map%volume
1593  else
1594  xi(1:params%dims(1)) = rdn(1:params%dims(1))
1595  x = map%get_x(xi)
1596  end if
1597  ! Maxwellian distribution of the temperature
1598  do i_v = 1,params%dims(2)
1599  call sll_s_normal_cdf_inv( rdn(i_v+params%dims(1)), 0.0_f64, 1.0_f64, &
1600  v(i_v))
1601  end do
1602  ! For multiple Gaussian, draw which one to take
1603  rnd_no = rdn(params%dims(1)+params%dims(2)+1)
1604  i_gauss = 1
1605  do while( rnd_no > delta(i_gauss) )
1606  i_gauss = i_gauss+1
1607  end do
1608  if( self%xiprofile ) then
1609  select type(p=>params)
1610  type is(sll_t_params_cos_gaussian_screwpinch)
1611  if( particle_group%species%q > 0._f64) then
1612  v(1:p%dims(2)) = v(1:p%dims(2)) * sqrt(p%profile%T_i(xi(1))/particle_group%species%m) + p%v_mean(:,i_gauss)
1613  else
1614  v(1:p%dims(2)) = v(1:p%dims(2)) * sqrt(p%profile%T_e(xi(1))/particle_group%species%m) + p%v_mean(:,i_gauss)
1615  end if
1616  class default
1617  ! Set weight according to value of perturbation
1618  wi(1) = params%eval_x_density(xi)*map%jacobian(xi)
1619  end select
1620  else
1621  ! Set weight according to value of perturbation
1622  wi(1) = params%eval_x_density(x)*map%jacobian(xi)
1623  v(1:params%dims(2)) = v(1:params%dims(2)) * params%v_thermal(:,i_gauss) + params%v_mean(:,i_gauss)
1624  end if
1625  elseif ( modulo(i_part, 2) == 0 ) then
1626  v(3) = - v(3) + 2.0_f64*params%v_mean(3,i_gauss)
1627  elseif ( modulo(i_part, 4) == 3 ) then
1628  v(2) = - v(2) + 2.0_f64*params%v_mean(2,i_gauss)
1629  v(1) = - v(1) + 2.0_f64*params%v_mean(1,i_gauss)
1630  end if
1631  ! Copy the generated numbers to the particle
1632  call particle_group%set_x(i_part, xi)
1633  call particle_group%set_v(i_part, v)
1634  ! Set weights.
1635  call particle_group%set_weights(i_part, &
1636  wi)
1637  ! w_local=w_local+wi(1)
1638  end do
1639 
1640  if( particle_group%species%q < 0._f64) then
1641  if( self%delta_perturb ) then
1642  call delta_function_perturbation( self, particle_group, [0._f64,0._f64,0._f64], [1._f64,1._f64,1._f64] )
1643  end if
1644  end if
1645 
1646  end subroutine sample_particle_sampling_sym_3d3v_trafo
1647 
1648  subroutine sll_s_particle_sampling_randomized_weights(particle_group, alpha, seed)
1649  class(sll_c_particle_group_base), target, intent(inout) :: particle_group
1650  sll_real64, intent(in) :: alpha
1651  sll_int32, intent(in), optional :: seed(:)
1652 
1653  sll_int32 :: i_part
1654  sll_real64 :: wi(particle_group%N_weights)
1655  sll_real64 :: pert, rnd(1)
1656 
1657  if (present(seed)) call random_seed(put=seed)
1658 
1659 
1660  do i_part = 1, particle_group%n_particles
1661  call random_number( rnd )
1662  wi = particle_group%get_weights(i_part)
1663  call sll_s_normal_cdf_inv( rnd(1), 0.0_f64, 1.0_f64, &
1664  pert)
1665  if (particle_group%n_weights > 2) then
1666  wi(3) = wi(3)-wi(1)
1667  end if
1668  wi(1) = wi(1) + alpha * pert
1669  if (particle_group%n_weights > 2) then
1670  wi(3) = wi(3)+wi(1)
1671  end if
1672  call particle_group%set_weights(i_part, wi)
1673  end do
1674 
1675 
1676  end subroutine sll_s_particle_sampling_randomized_weights
1677 
1678 
1679 end module sll_m_particle_sampling
sll_m_collective::sll_o_collective_allreduce
Combines values from all processes and distributes the result back to all processes.
Definition: sll_m_collective.F90:312
sll_m_collective
Parallelizing facility.
Definition: sll_m_collective.F90:128
sll_m_collective::sll_v_world_collective
type(sll_t_collective_t), pointer, public sll_v_world_collective
Control of subset assignment. Processes with the same color are in the same new communicator.
Definition: sll_m_collective.F90:255
sll_m_control_variate
Control variate.
Definition: sll_m_control_variate.F90:6
sll_m_initial_distribution
Parameters to define common initial distributions.
Definition: sll_m_initial_distribution.F90:5
sll_m_mapping_3d
Module interfaces for coordinate transformation.
Definition: sll_m_mapping_3d.F90:13
sll_m_particle_group_base
Definition: sll_m_particle_group_base.F90:1
sll_m_particle_sampling
Particle initializer class with various functions to initialize a particle.
Definition: sll_m_particle_sampling.F90:5
sll_m_particle_sampling::delta_function_perturbation
subroutine delta_function_perturbation(self, particle_group, xmin, Lx)
Definition: sll_m_particle_sampling.F90:1388
sll_m_particle_sampling::sample_particle_sampling_sym_1d2v
subroutine sample_particle_sampling_sym_1d2v(self, particle_group, params, xmin, Lx)
Helper function for antithetic sampling in 1d2v.
Definition: sll_m_particle_sampling.F90:616
sll_m_particle_sampling::sample_particle_sampling_sym_uni_3d3v
subroutine sample_particle_sampling_sym_uni_3d3v(self, particle_group, params, xmin, Lx)
Helper function for antithetic sampling in 1d2v.
Definition: sll_m_particle_sampling.F90:1159
sll_m_particle_sampling::sll_p_symmetric_all
integer(kind=i32), parameter sll_p_symmetric_all
2^(dim_x+dim_v) points constructed from each drawn particle, all possible combinations of reflections...
Definition: sll_m_particle_sampling.F90:54
sll_m_particle_sampling::sample_particle_sampling_unix
subroutine sample_particle_sampling_unix(self, particle_group, params, xmin, Lx)
Helper function for antithetic sampling in 1d2v.
Definition: sll_m_particle_sampling.F90:804
sll_m_particle_sampling::sample_particle_sampling_sym_1d2v_trafo
subroutine sample_particle_sampling_sym_1d2v_trafo(self, particle_group, params, Lx, map)
Helper function for antithetic sampling in 1d2v.
Definition: sll_m_particle_sampling.F90:1419
sll_m_particle_sampling::sample_cv_particle_sampling
subroutine sample_cv_particle_sampling(self, particle_group, params, xmin, Lx, control_variate, time, map, lindf)
Sample with control variate (we assume that the particle weights are given in the following order: (f...
Definition: sll_m_particle_sampling.F90:259
sll_m_particle_sampling::init_particle_sampling
subroutine init_particle_sampling(self, sampling_type, dims, n_particles_local, rank, delta_perturb, delta_eps)
Initializer.
Definition: sll_m_particle_sampling.F90:96
sll_m_particle_sampling::sll_p_random_numbers
integer(kind=i32), parameter sll_p_random_numbers
draw random numbers
Definition: sll_m_particle_sampling.F90:50
sll_m_particle_sampling::sll_p_uniformx_negative
integer(kind=i32), parameter sll_p_uniformx_negative
Definition: sll_m_particle_sampling.F90:56
sll_m_particle_sampling::sample_particle_sampling_all
subroutine sample_particle_sampling_all(self, particle_group, params, xmin, Lx)
Helper function for pure sampling.
Definition: sll_m_particle_sampling.F90:504
sll_m_particle_sampling::sll_p_symmetric_negative
integer(kind=i32), parameter sll_p_symmetric_negative
one additional point for each drawd point, the one reflected along all directions
Definition: sll_m_particle_sampling.F90:55
sll_m_particle_sampling::sll_p_standard
integer(kind=i32), parameter sll_p_standard
each particle is drawn separately
Definition: sll_m_particle_sampling.F90:53
sll_m_particle_sampling::sample_particle_sampling_all_trafo
subroutine sample_particle_sampling_all_trafo(self, particle_group, params, xmin, Lx, map)
Helper function for pure sampling.
Definition: sll_m_particle_sampling.F90:1254
sll_m_particle_sampling::sample_particle_sampling_unix_1d2v
subroutine sample_particle_sampling_unix_1d2v(self, particle_group, params, xmin, Lx)
Helper function for antithetic sampling in 1d2v.
Definition: sll_m_particle_sampling.F90:939
sll_m_particle_sampling::sample_particle_sampling
subroutine sample_particle_sampling(self, particle_group, params, xmin, Lx, map)
Sample from distribution defined by params.
Definition: sll_m_particle_sampling.F90:444
sll_m_particle_sampling::sample_particle_sampling_sym_3d3v_trafo
subroutine sample_particle_sampling_sym_3d3v_trafo(self, particle_group, params, xmin, Lx, map)
Helper function for antithetic sampling in 1d2v.
Definition: sll_m_particle_sampling.F90:1516
sll_m_particle_sampling::sample_particle_sampling_sym_3d3v
subroutine sample_particle_sampling_sym_3d3v(self, particle_group, params, xmin, Lx)
Helper function for antithetic sampling in 3d3v.
Definition: sll_m_particle_sampling.F90:1053
sll_m_particle_sampling::sample_particle_sampling_sym2_1d2v
subroutine sample_particle_sampling_sym2_1d2v(self, particle_group, params, xmin, Lx)
Helper function for antithetic sampling in 1d2v.
Definition: sll_m_particle_sampling.F90:714
sll_m_particle_sampling::reset_seed_jump
subroutine reset_seed_jump(self, jump)
Definition: sll_m_particle_sampling.F90:243
sll_m_particle_sampling::sll_s_particle_sampling_randomized_weights
subroutine, public sll_s_particle_sampling_randomized_weights(particle_group, alpha, seed)
Definition: sll_m_particle_sampling.F90:1649
sll_m_particle_sampling::sll_p_sobol_numbers
integer(kind=i32), parameter sll_p_sobol_numbers
draw sobol numbers
Definition: sll_m_particle_sampling.F90:51
sll_m_particle_sampling::free_particle_sampling
subroutine free_particle_sampling(self)
Descructor.
Definition: sll_m_particle_sampling.F90:88
sll_m_control_variate::sll_t_control_variate
Control variate object.
Definition: sll_m_control_variate.F90:23
sll_m_initial_distribution::sll_c_distribution_params
Abstract data type for parameters of initial distribution.
Definition: sll_m_initial_distribution.F90:48
sll_m_initial_distribution::sll_t_params_cos_gaussian_screwpinch
Definition: sll_m_initial_distribution.F90:100
sll_m_initial_distribution::sll_t_params_cos_gaussian
Data type for distribution function with (multiple) Gaussians in v and one plus cosine perturbations ...
Definition: sll_m_initial_distribution.F90:65
sll_m_initial_distribution::sll_t_params_noise_gaussian
Definition: sll_m_initial_distribution.F90:81
sll_m_mapping_3d::sll_t_mapping_3d
type collecting functions for analytical coordinate mapping
Definition: sll_m_mapping_3d.F90:68
sll_m_particle_group_base::sll_c_particle_group_base
Definition: sll_m_particle_group_base.F90:42
sll_m_particle_sampling::sll_t_particle_sampling
Data type for particle sampling.
Definition: sll_m_particle_sampling.F90:59
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:40 for Semi-Lagrangian Library by   doxygen 1.9.1