Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_particle_mesh_coupling_spline_strong_1d.F90
Go to the documentation of this file.
1 
6 module sll_m_particle_mesh_coupling_spline_strong_1d
7 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8 #include "sll_assert.h"
9 #include "sll_memory.h"
10 #include "sll_working_precision.h"
11 #include "sll_errors.h"
12 
13  use sll_m_low_level_bsplines, only: &
14  sll_s_uniform_bsplines_eval_basis
15 
16  use sll_m_gauss_legendre_integration, only : &
17  sll_f_gauss_legendre_points_and_weights
18 
19  use sll_m_particle_mesh_coupling_base_1d, only: &
20  sll_p_collocation, &
21  sll_c_particle_mesh_coupling_1d, &
22  sll_p_galerkin
23 
24  use sll_m_particle_group_base, only: &
25  sll_c_particle_group_base
26 
27  use sll_m_splines_pp
28 
29  implicit none
30 
31  public :: &
32  sll_t_particle_mesh_coupling_spline_strong_1d, &
33  sll_s_new_particle_mesh_coupling_spline_strong_1d
34 
35  private
36 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
38  type, extends(sll_c_particle_mesh_coupling_1d) :: sll_t_particle_mesh_coupling_spline_strong_1d
39  ! Member variables form abstract type
40  logical :: integ
41  logical :: eval_same = .true.
42  sll_int32 :: index_shift
43 
44  ! Information about the 1d mesh
45  sll_int32 :: n_span
46  type(sll_t_spline_pp_1d) :: spline_pp
47 
48 
49  sll_real64, allocatable :: spline_val(:)
50  sll_real64, allocatable :: spline_val_more(:)
51 
52  sll_real64, allocatable :: spline_val_prim(:)
53  sll_real64, allocatable :: spline_val_prim_old(:)
54  sll_real64, allocatable :: spline_val_prim_new(:)
55 
56  sll_int32 :: n_quad_points
57  sll_real64, allocatable :: quad_xw(:,:)
58 
59  contains
60 
61  procedure :: add_charge => add_charge_single_spline_strong_1d
62  procedure :: add_particle_mass => add_particle_mass_spline_strong_1d
63  procedure :: add_particle_mass_full => add_particle_mass_spline_strong_1d
64  procedure :: evaluate => evaluate_field_single_spline_strong_1d
65  procedure :: evaluate_multiple => evaluate_multiple_spline_strong_1d
66  procedure :: add_current_update_v => add_current_update_v_spline_strong_1d_quadrature
67  !procedure :: add_current_update_v_pp => add_current_update_v_pp_spline_strong_1d !> Add contribution of pne particle to the current density and update velocity
68  procedure :: add_current => add_current_spline_strong_1d
69  procedure :: add_current_evaluate => add_current_evaluate_spline_strong_1d_quadrature
70 
71  procedure :: init => init_spline_strong_1d
72  procedure :: free => free_spline_strong_1d
73 
74 
75  end type sll_t_particle_mesh_coupling_spline_strong_1d
76 
77 contains
78 
79  ! TODOTODOTODOTODO
80 
81  ! Check the /delta_x scaling in the add_charge function
82  subroutine helper_normalized_position ( self, position, index, axi)
83  class( sll_t_particle_mesh_coupling_spline_strong_1d ), intent(inout) :: self
84  sll_real64, intent( in ) :: position(self%dim)
85  sll_int32, intent(out) :: index
86  sll_real64, intent(out) :: axi(self%dim)
87 
88  sll_real64 :: xi(1)
89 
90  xi(1) = (position(1) - self%domain(1))/self%delta_x
91  index = floor(xi(1))+1
92  xi(1) = xi(1) - real(index-1, f64)
93 
94  ! Old version before changed to option of evaluating in other points than the grid point
95 !!$ if ( modulo(self%spline_degree, 2) == 0 ) then
96 !!$ if (xi(1) < 0.5_f64 ) then
97 !!$ axi(1) = 0.5_f64 - xi(1)
98 !!$ index = index - self%spline_degree/2
99 !!$ else
100 !!$ axi(1) = 1.5_f64 - xi(1)
101 !!$ index = index - self%spline_degree/2 + 1
102 !!$ end if
103 !!$ else
104 !!$ axi(1) = 1.0_f64 - xi(1)
105 !!$ index = index - (self%n_span)/2+1
106 !!$ end if
107  if ( self%eval_same .eqv. .true. ) then ! case odd spline with evaluation at grid points or even spline with evaluation at mid points
108  axi(1) = 1.0_f64 - xi(1)
109  index = index - self%index_shift!(self%spline_degree)/2
110  else ! case even spline with evaluation at grid points or odd spline with evaluation at mid points
111  if (xi(1) < 0.5_f64 ) then
112  axi(1) = 0.5_f64 - xi(1)
113  index = index - self%index_shift!(self%n_span-1)/2
114  else
115  axi(1) = 1.5_f64 - xi(1)
116  index = index - self%index_shift + 1!(self%n_span-1)/2 + 1
117  end if
118  end if
119 
120  ! print*, axi, index
121  ! print*, (self%spline_degree-1)/2, (self%spline_degree)/2
122 
123  end subroutine helper_normalized_position
124 
125 
126  !---------------------------------------------------------------------------!
128  subroutine add_charge_single_spline_strong_1d(self, position, marker_charge, rho_dofs)
129  class( sll_t_particle_mesh_coupling_spline_strong_1d ), intent(inout) :: self
130  sll_real64, intent( in ) :: position(self%dim)
131  sll_real64, intent( in ) :: marker_charge
132  sll_real64, intent( inout ) :: rho_dofs(self%n_dofs)
133 
134  ! NOTE THAT THIS IS THE VARIANT FOR INT=TRUE ONLY
135  !local variables
136  sll_int32 :: i1, i
137  sll_int32 :: index1d, index
138  sll_real64 :: axi(1)
139 
140  call helper_normalized_position( self, position, index, axi )
141 
142  if ( self%integ .eqv. .false. ) then
143  ! Evaluate the spline function
144  call sll_s_uniform_bsplines_eval_basis( self%spline_degree, axi(1), &
145  self%spline_val )
146 
147  do i1 = 1, self%n_span
148  index1d = modulo(index+i1-2, self%n_cells)+1
149  rho_dofs(index1d) = rho_dofs(index1d) + &
150  marker_charge * self%spline_val(self%n_span-i1+1)/ self%delta_x
151  end do
152  else
153 
154  ! Old version with primitive formula
155  ! Evaluate the primitive function
156  do i=1, self%n_span
157  self%spline_val_prim(self%n_span-i+2) = sll_f_spline_pp_horner_1d(self%spline_degree+1, self%spline_pp%poly_coeffs_fpa, axi(1), i)
158  end do
159 
160  ! fill the right-hand-side
161  do i=1, self%n_span+1
162  index1d = modulo(index+i-2,self%n_cells)+1
163  rho_dofs(index1d) = rho_dofs(index1d) + &
164  marker_charge * ( self%spline_val_prim(i+1) - self%spline_val_prim(i) )/ self%delta_x
165  end do
166  end if
167 
168  end subroutine add_charge_single_spline_strong_1d
169 
171  subroutine add_current_update_v_spline_strong_1d_quadrature (self, position_old, position_new, marker_charge, qoverm, bfield_dofs, vi, j_dofs)
172  class(sll_t_particle_mesh_coupling_spline_strong_1d), intent(inout) :: self
173  sll_real64, intent(in) :: position_old(self%dim)
174  sll_real64, intent(in) :: position_new(self%dim)
175  sll_real64, intent(in) :: marker_charge
176  sll_real64, intent(in) :: qoverm
177  sll_real64, intent(in) :: bfield_dofs(:)
178  sll_real64, intent(inout) :: vi(:)
179  sll_real64, intent(inout) :: j_dofs(:)
180 
181  !local variables
182  sll_int32 :: ind
183  sll_int32 :: index_new, index_old
184  sll_real64 :: axi_new(1), axi_old(1)
185 
186  call helper_normalized_position( self, position_old, index_old, axi_old )
187  call helper_normalized_position( self, position_new, index_new, axi_new )
188 
189  if (index_old == index_new) then
190  call current_v_local( self, axi_old(1), axi_new(1), index_old, marker_charge, qoverm, vi(2), j_dofs, bfield_dofs )
191  elseif(index_old < index_new ) then
192  call current_v_local( self, axi_old(1), 0.0_f64, index_old, marker_charge, qoverm, vi(2), j_dofs, bfield_dofs )
193  do ind = index_old+1, index_new-1
194  call current_v_local( self, 1.0_f64, 0.0_f64, ind , marker_charge, qoverm, vi(2), j_dofs, bfield_dofs)
195  end do
196  call current_v_local( self, 1.0_f64, axi_new(1), index_new, marker_charge, qoverm, vi(2), j_dofs, bfield_dofs )
197  else
198  call current_v_local( self, axi_old(1), 1.0_f64, index_old, marker_charge, qoverm, vi(2), j_dofs, bfield_dofs )
199  do ind = index_new+1, index_old-1
200  call current_v_local( self, 0.0_f64, 1.0_f64, ind , marker_charge, qoverm, vi(2), j_dofs, bfield_dofs)
201  end do
202  call current_v_local( self, 0.0_f64, axi_new(1), index_new, marker_charge, qoverm, vi(2), j_dofs, bfield_dofs )
203  end if
204 
205 
206  end subroutine add_current_update_v_spline_strong_1d_quadrature
207 
208 
209  subroutine current_v_local( self, upper, lower, box, marker_charge, qoverm, vi, j_dofs, bfield_dofs)
210  class(sll_t_particle_mesh_coupling_spline_strong_1d), intent(inout) :: self
211  sll_real64, intent(in) :: lower
212  sll_real64, intent(in) :: upper
213  sll_int32, intent(in) :: box
214  sll_real64, intent(in) :: marker_charge
215  sll_real64, intent(in) :: qoverm
216  sll_real64, intent(inout) :: vi
217  sll_real64, intent(in) :: bfield_dofs(self%n_dofs)
218  sll_real64, intent(inout) :: j_dofs(self%n_dofs)
219 
220  !Local variables
221  sll_int32 :: ind, i_grid, i_mod, n_cells, j
222  sll_real64 :: c1, c2
223  n_cells = self%n_cells
224 
225  c1 = 0.5_f64*(upper-lower)
226  c2 = 0.5_f64*(upper+lower)
227 
228  call sll_s_uniform_bsplines_eval_basis(self%spline_degree, c1*self%quad_xw(1,1)+c2, &
229  self%spline_val)
230  self%spline_val = self%spline_val * (self%quad_xw(2,1)*c1)
231  do j=2,self%n_quad_points
232  call sll_s_uniform_bsplines_eval_basis(self%spline_degree, c1*self%quad_xw(1,j)+c2, &
233  self%spline_val_more)
234  self%spline_val = self%spline_val + self%spline_val_more * (self%quad_xw(2,j)*c1)
235  end do
236 ! self%spline_val = self%spline_val * (self%delta_x)
237 ! do i=1, self%n_span
238 ! self%spline_val_prim_old(self%n_span-i+1) = sll_f_spline_pp_horner_1d(self%spline_degree+1, self%spline_pp%poly_coeffs_fpa, upper, i)
239 ! self%spline_val_prim_new(self%n_span-i+1) = sll_f_spline_pp_horner_1d(self%spline_degree+1, self%spline_pp%poly_coeffs_fpa, lower, i)
240 ! end do
241 
242  ind = self%spline_degree+1
243  do i_grid = box, box+self%spline_degree
244  i_mod = modulo(i_grid-1, n_cells ) +1
245  j_dofs(i_mod) = j_dofs(i_mod) + marker_charge*self%spline_val(ind)
246  vi = vi - qoverm * self%spline_val(ind)*bfield_dofs(i_mod)
247  ind = ind-1
248  end do
249 
250 
251  end subroutine current_v_local
252 
253 
255  subroutine add_current_update_v_spline_strong_1d (self, position_old, position_new, marker_charge, qoverm, bfield_dofs, vi, j_dofs)
256  class(sll_t_particle_mesh_coupling_spline_strong_1d), intent(inout) :: self
257  sll_real64, intent(in) :: position_old(self%dim)
258  sll_real64, intent(in) :: position_new(self%dim)
259  sll_real64, intent(in) :: marker_charge
260  sll_real64, intent(in) :: qoverm
261  sll_real64, intent(in) :: bfield_dofs(:)
262  sll_real64, intent(inout) :: vi(:)
263  sll_real64, intent(inout) :: j_dofs(:)
264  !local variables
265  sll_int32 :: i, i1, diff
266  sll_int32 :: index_new, index_old
267  sll_real64 :: axi_new(1), axi_old(1)
268 
269  call helper_normalized_position( self, position_old, index_old, axi_old )
270  call helper_normalized_position( self, position_new, index_new, axi_new )
271 
272  ! Evaluate the primitive function
273  do i=1, self%n_span
274  self%spline_val_prim_old(self%n_span-i+1) = sll_f_spline_pp_horner_1d(self%spline_degree+1, self%spline_pp%poly_coeffs_fpa, axi_old(1), i)
275  self%spline_val_prim_new(self%n_span-i+1) = sll_f_spline_pp_horner_1d(self%spline_degree+1, self%spline_pp%poly_coeffs_fpa, axi_new(1), i)
276  end do
277 
278  !call sll_s_spline_pp_horner_primitive_1d(self%spline_val, self%spline_degree, self%spline_pp%poly_coeffs_fp, lower)
279  !call sll_s_spline_pp_horner_primitive_1d(self%spline_val, self%spline_degree, self%spline_pp%poly_coeffs_fp, lower)
280  !print*, axi_old(1), axi_new(1)
281  !print*, self%spline_val_prim_new
282  !print*, self%spline_val_prim_old
283  !stop
284  ! do i=1, self%n_span
285  ! self%spline_val_prim_old(i) = sll_f_spline_pp_horner_1d(self%spline_degree, self%spline_pp%poly_coeffs, axi_old(1), i) / self%delta_x * (position_new(1) - position_old(1))
286  ! self%spline_val_prim_new(i) = -sll_f_spline_pp_horner_1d(self%spline_degree, self%spline_pp%poly_coeffs, axi_new(1), i) / self%delta_x * (position_new(1)-position_old(1))
287  ! end do
288 
289  if (index_old<index_new) then
290  diff = index_new-index_old
291 
292  do i=0,min(diff-1, self%spline_degree)
293  i1 = modulo(index_old+i-1, self%n_cells) + 1
294  j_dofs(i1) = j_dofs(i1) + self%spline_val_prim_old(i+1)* marker_charge
295  vi(2) = vi(2) - qoverm * bfield_dofs(i1) * ( self%spline_val_prim_old(i+1))
296  end do
297  do i = self%spline_degree+1, diff-1
298  i1 = modulo(index_old+i-1, self%n_cells) + 1
299  j_dofs(i1) = j_dofs(i1) + marker_charge
300  vi(2) = vi(2) - qoverm * bfield_dofs(i1)
301  end do
302  do i = diff, self%spline_degree
303  i1 = modulo(index_old+i-1, self%n_cells) + 1
304  j_dofs(i1) = j_dofs(i1) + ( self%spline_val_prim_old(i+1) - self%spline_val_prim_new(i+1-diff)) &
305  * marker_charge
306  vi(2) = vi(2) - qoverm * bfield_dofs(i1) * ( self%spline_val_prim_old(i+1) - self%spline_val_prim_new(i+1-diff))
307  end do
308  do i=max(self%spline_degree+1, diff), self%spline_degree+diff
309  i1 = modulo(index_old+i-1, self%n_cells) + 1
310  j_dofs(i1) = j_dofs(i1) + ( 1.0_f64 - self%spline_val_prim_new(i+1-diff))*marker_charge
311  vi(2) = vi(2) - qoverm * bfield_dofs(i1) * ( 1.0_f64 - self%spline_val_prim_new(i+1-diff))
312  end do
313  else
314  diff = index_old-index_new
315  do i=0,min(diff-1, self%spline_degree)
316  i1 = modulo(index_new+i-1, self%n_cells) + 1
317  j_dofs(i1) = j_dofs(i1) - self%spline_val_prim_new(i+1)*marker_charge
318  vi(2) = vi(2) - qoverm * bfield_dofs(i1) * ( - self%spline_val_prim_new(i+1))
319  end do
320  do i=self%spline_degree+1, diff-1
321  i1 = modulo(index_new+i-1, self%n_cells) + 1
322  j_dofs(i1) = j_dofs(i1) - marker_charge
323  vi(2) = vi(2) + qoverm * bfield_dofs(i1)
324  end do
325 
326  do i = diff, self%spline_degree
327  i1 = modulo(index_new+i-1, self%n_cells) + 1
328  j_dofs(i1) = j_dofs(i1) +( self%spline_val_prim_old(i+1-diff) - self%spline_val_prim_new(i+1)) &
329  * marker_charge
330  vi(2) = vi(2) - qoverm * bfield_dofs(i1) * ( self%spline_val_prim_old(i+1-diff) - self%spline_val_prim_new(i+1))
331  end do
332  do i=max(self%spline_degree+1, diff), self%spline_degree+diff
333  i1 = modulo(index_new+i-1, self%n_cells) + 1
334  j_dofs(i1) = j_dofs(i1) + (-1.0_f64 + self%spline_val_prim_old(i+1-diff)) &
335  * marker_charge
336  vi(2) = vi(2) - qoverm * bfield_dofs(i1) * ( self%spline_val_prim_old(i+1-diff) - 1.0_f64 )
337  end do
338  end if
339 
340  end subroutine add_current_update_v_spline_strong_1d
341 
342 
343  subroutine evaluate_field_single_spline_strong_1d(self, position, field_dofs, field_value)
344  class(sll_t_particle_mesh_coupling_spline_strong_1d), intent( inout ) :: self
345  sll_real64, intent( in ) :: position(self%dim)
346  sll_real64, intent( in ) :: field_dofs(self%n_dofs)
347  sll_real64, intent( out ) :: field_value
348 
349  !local variables
350  sll_int32 :: i1, i
351  sll_int32 :: index1d, index
352  sll_real64 :: axi(1)
353 
354 
355  call helper_normalized_position( self, position, index, axi )
356 
357  if ( self%integ .eqv. .true. ) then
358  ! Evaluate the primitive function
359  do i=1, self%n_span
360  self%spline_val_prim(self%n_span-i+2) = sll_f_spline_pp_horner_1d(self%spline_degree+1, self%spline_pp%poly_coeffs_fpa, axi(1), i)
361  end do
362 
363  field_value = 0.0_f64
364  do i1 = 1, self%n_span+1
365  index1d = modulo(index+i1-2, self%n_cells)+1
366  field_value = field_value + &
367  field_dofs(index1d) * &
368  (self%spline_val_prim(i1+1)-self%spline_val_prim(i1)) / &
369  self%delta_x
370  end do
371  else
372  ! Evaluate the spline function
373  call sll_s_uniform_bsplines_eval_basis( self%spline_degree, axi(1), &
374  self%spline_val )
375 
376  field_value = 0.0_f64
377  do i1 = 1, self%n_span
378  index1d = modulo(index+i1-2, self%n_cells)+1
379  field_value = field_value + &
380  field_dofs(index1d) * &
381  self%spline_val(self%n_span-i1+1) / &
382  self%delta_x
383  end do
384 
385 
386  end if
387 
388  end subroutine evaluate_field_single_spline_strong_1d
389 
390 
392  subroutine add_current_evaluate_spline_strong_1d_quadrature (self, position_old, position_new, marker_charge, vbar, field_dofs, j_dofs, field)
393  class(sll_t_particle_mesh_coupling_spline_strong_1d), intent(inout) :: self
394  sll_real64, intent(in) :: position_old(self%dim)
395  sll_real64, intent(in) :: position_new(self%dim)
396  sll_real64, intent(in) :: marker_charge
397  sll_real64, intent( in ) :: vbar
398  sll_real64, intent(in) :: field_dofs(self%n_dofs)
399  sll_real64, intent(inout) :: j_dofs(self%n_dofs)
400  sll_real64, intent( out ) :: field
401 
402  !local variables
403  sll_int32 :: ind
404  sll_int32 :: index_new, index_old
405  sll_real64 :: axi_new(1), axi_old(1)
406 
407 
408  field = 0.0_f64
409  call helper_normalized_position( self, position_old, index_old, axi_old )
410  call helper_normalized_position( self, position_new, index_new, axi_new )
411 
412  if (index_old == index_new) then
413  call current_eval_local( self, axi_old(1), axi_new(1), index_old, marker_charge, field, j_dofs, field_dofs )
414  elseif(index_old < index_new ) then
415  call current_eval_local( self, axi_old(1), 0.0_f64, index_old, marker_charge, field, j_dofs, field_dofs )
416  do ind = index_old+1, index_new-1
417  call current_eval_local( self, 1.0_f64, 0.0_f64, ind , marker_charge, field, j_dofs, field_dofs)
418  end do
419  call current_eval_local( self, 1.0_f64, axi_new(1), index_new, marker_charge, field, j_dofs, field_dofs )
420  else
421  call current_eval_local( self, axi_old(1), 1.0_f64, index_old, marker_charge, field, j_dofs, field_dofs )
422  do ind = index_new+1, index_old-1
423  call current_eval_local( self, 0.0_f64, 1.0_f64, ind , marker_charge, field, j_dofs, field_dofs)
424  end do
425  call current_eval_local( self, 0.0_f64, axi_new(1), index_new, marker_charge, field, j_dofs, field_dofs )
426  end if
427 
428  field = field/vbar
429 
430  end subroutine add_current_evaluate_spline_strong_1d_quadrature
431 
432 
433 
434  subroutine current_eval_local( self, upper, lower, box, marker_charge, field, j_dofs, field_dofs)
435  class(sll_t_particle_mesh_coupling_spline_strong_1d), intent(inout) :: self
436  sll_real64, intent(in) :: lower
437  sll_real64, intent(in) :: upper
438  sll_int32, intent(in) :: box
439  sll_real64, intent(in) :: marker_charge
440  sll_real64, intent(inout) :: field
441  sll_real64, intent(in) :: field_dofs(self%n_dofs)
442  sll_real64, intent(inout) :: j_dofs(self%n_dofs)
443 
444  !Local variables
445  sll_int32 :: ind, i_grid, i_mod, n_cells, j
446  sll_real64 :: c1, c2
447  n_cells = self%n_cells
448 
449  c1 = 0.5_f64*(upper-lower)
450  c2 = 0.5_f64*(upper+lower)
451 
452  call sll_s_uniform_bsplines_eval_basis(self%spline_degree, c1*self%quad_xw(1,1)+c2, &
453  self%spline_val)
454  self%spline_val = self%spline_val * (self%quad_xw(2,1)*c1)
455  do j=2,self%n_quad_points
456  call sll_s_uniform_bsplines_eval_basis(self%spline_degree, c1*self%quad_xw(1,j)+c2, &
457  self%spline_val_more)
458  self%spline_val = self%spline_val + self%spline_val_more * (self%quad_xw(2,j)*c1)
459  end do
460 ! self%spline_val = self%spline_val * (self%delta_x)
461 ! do i=1, self%n_span
462 ! self%spline_val_prim_old(self%n_span-i+1) = sll_f_spline_pp_horner_1d(self%spline_degree+1, self%spline_pp%poly_coeffs_fpa, upper, i)
463 ! self%spline_val_prim_new(self%n_span-i+1) = sll_f_spline_pp_horner_1d(self%spline_degree+1, self%spline_pp%poly_coeffs_fpa, lower, i)
464 ! end do
465 
466  ind = self%spline_degree+1
467  do i_grid = box, box+self%spline_degree
468  i_mod = modulo(i_grid-1, n_cells ) +1
469  j_dofs(i_mod) = j_dofs(i_mod) + marker_charge*self%spline_val(ind)
470  field = field + field_dofs(i_mod)* self%spline_val(ind)
471  ind = ind-1
472  end do
473 
474 
475  end subroutine current_eval_local
476 
477 
478 
479 
480  subroutine add_current_evaluate_spline_strong_1d(self, position_old, position_new, marker_charge, vbar, field_dofs, j_dofs, field)
481  class( sll_t_particle_mesh_coupling_spline_strong_1d ), intent(inout) :: self
482  sll_real64, intent( in ) :: position_old(self%dim)
483  sll_real64, intent( in ) :: position_new(self%dim)
484  sll_real64, intent( in ) :: marker_charge
485  sll_real64, intent( in ) :: vbar
486  sll_real64, intent( in ) :: field_dofs(self%n_dofs)
487  sll_real64, intent( inout ) :: j_dofs(self%n_dofs)
488  sll_real64, intent( out ) :: field
489  !local variables
490  sll_int32 :: i, i1, diff
491  sll_int32 :: index_new, index_old
492  sll_real64 :: axi_new(1), axi_old(1)
493 
494  call helper_normalized_position( self, position_old, index_old, axi_old )
495  call helper_normalized_position( self, position_new, index_new, axi_new )
496 
497  ! Evaluate the primitive function
498  do i=1, self%n_span
499  self%spline_val_prim_old(self%n_span-i+1) = sll_f_spline_pp_horner_1d(self%spline_degree+1, self%spline_pp%poly_coeffs_fpa, axi_old(1), i)
500  self%spline_val_prim_new(self%n_span-i+1) = sll_f_spline_pp_horner_1d(self%spline_degree+1, self%spline_pp%poly_coeffs_fpa, axi_new(1), i)
501  end do
502 
503  field = 0.0_f64
504 
505  if (index_old<index_new) then
506  diff = index_new-index_old
507 
508  do i=0,min(diff-1, self%spline_degree)
509  i1 = modulo(index_old+i-1, self%n_cells) + 1
510  j_dofs(i1) = j_dofs(i1) + self%spline_val_prim_old(i+1)* marker_charge
511  field = field + field_dofs(i1) * ( self%spline_val_prim_old(i+1))
512  end do
513  do i = self%spline_degree+1, diff-1
514  i1 = modulo(index_old+i-1, self%n_cells) + 1
515  j_dofs(i1) = j_dofs(i1) + marker_charge
516  field = field + field_dofs(i1)
517  end do
518  do i = diff, self%spline_degree
519  i1 = modulo(index_old+i-1, self%n_cells) + 1
520  j_dofs(i1) = j_dofs(i1) + ( self%spline_val_prim_old(i+1) - self%spline_val_prim_new(i+1-diff)) &
521  * marker_charge
522  field = field + field_dofs(i1) * ( self%spline_val_prim_old(i+1) - self%spline_val_prim_new(i+1-diff))
523  end do
524  do i=max(self%spline_degree+1, diff), self%spline_degree+diff
525  i1 = modulo(index_old+i-1, self%n_cells) + 1
526  j_dofs(i1) = j_dofs(i1) + ( 1.0_f64 - self%spline_val_prim_new(i+1-diff))*marker_charge
527  field = field + field_dofs(i1) * ( 1.0_f64 - self%spline_val_prim_new(i+1-diff))
528  end do
529  else
530  diff = index_old-index_new
531  do i=0,min(diff-1, self%spline_degree)
532  i1 = modulo(index_new+i-1, self%n_cells) + 1
533  j_dofs(i1) = j_dofs(i1) - self%spline_val_prim_new(i+1)*marker_charge
534  field = field + field_dofs(i1) * ( - self%spline_val_prim_new(i+1))
535  end do
536  do i=self%spline_degree+1, diff-1
537  i1 = modulo(index_new+i-1, self%n_cells) + 1
538  j_dofs(i1) = j_dofs(i1) - marker_charge
539  field = field - field_dofs(i1)
540  end do
541 
542  do i = diff, self%spline_degree
543  i1 = modulo(index_new+i-1, self%n_cells) + 1
544  j_dofs(i1) = j_dofs(i1) +( self%spline_val_prim_old(i+1-diff) - self%spline_val_prim_new(i+1)) &
545  * marker_charge
546  field = field + field_dofs(i1) * ( self%spline_val_prim_old(i+1-diff) - self%spline_val_prim_new(i+1))
547  end do
548  do i=max(self%spline_degree+1, diff), self%spline_degree+diff
549  i1 = modulo(index_new+i-1, self%n_cells) + 1
550  j_dofs(i1) = j_dofs(i1) + (-1.0_f64 + self%spline_val_prim_old(i+1-diff)) &
551  * marker_charge
552  field = field + field_dofs(i1) * ( self%spline_val_prim_old(i+1-diff) - 1.0_f64 )
553  end do
554  end if
555 
556  field = field/vbar
557 
558  end subroutine add_current_evaluate_spline_strong_1d
559 
560 
562  subroutine init_spline_strong_1d( self, domain, n_cells, spline_degree, integ, eval_grid_points )
563  class(sll_t_particle_mesh_coupling_spline_strong_1d), intent( inout ) :: self
564  sll_int32, intent(in) :: n_cells
565  sll_real64, intent(in) :: domain(2)
566  sll_int32, intent(in) :: spline_degree
567  logical, intent( in ) :: integ
568  logical, intent( in ) :: eval_grid_points
569 
570  !local variables
571  sll_int32 :: ierr
572 
573  self%integ = .false.!integ
574  self%dim = 1
575 
576  ! Store grid information
577  self%domain = domain
578  self%n_cells = n_cells
579  self%n_dofs = n_cells
580  self%delta_x = (self%domain(2)-self%domain(1))/real(n_cells, f64)
581 
582 
583  ! Initialize information on the spline
584  self%spline_degree = spline_degree
585  self%n_span = spline_degree + 1
586 
587 ! self%n_quad_points = (self%spline_degree+2)/2
588 
589  ALLOCATE( self%spline_val(self%n_span), stat = ierr)
590  sll_assert( ierr == 0 )
591  ALLOCATE( self%spline_val_more(self%n_span), stat = ierr)
592  sll_assert( ierr == 0 )
593  ALLOCATE( self%spline_val_prim(self%n_span+2), stat = ierr )
594  sll_assert( ierr == 0 )
595  self%spline_val_prim(1) = 0.0_f64
596  self%spline_val_prim(self%n_span+2) = 1.0_f64
597  ALLOCATE( self%spline_val_prim_old(self%n_span), stat = ierr )
598  sll_assert( ierr == 0 )
599  ALLOCATE( self%spline_val_prim_new(self%n_span), stat = ierr )
600  sll_assert( ierr == 0 )
601 
602  if ( self%spline_degree < 6 ) then
603  call sll_s_spline_pp_init_1d( self%spline_pp, spline_degree, n_cells)
604  end if
618  if ( integ .eqv. .false. ) then
619  if ( eval_grid_points .eqv. .true. ) then
620  if ( modulo( self%spline_degree, 2 ) == 0 ) then
621  self%eval_same = .false.
622  self%index_shift = (self%spline_degree)/2
623  else
624  self%eval_same = .true.
625  self%index_shift = (self%spline_degree-1)/2
626  end if
627  else
628  if ( modulo( self%spline_degree, 2 ) == 0 ) then
629  self%eval_same = .true.
630  self%index_shift = (self%spline_degree)/2
631  else
632  self%eval_same = .false.
633  self%index_shift = (self%spline_degree+1)/2
634  end if
635  end if
636  else
637  if ( eval_grid_points .eqv. .true. ) then
638  if ( modulo( self%spline_degree, 2 ) == 0 ) then
639  self%eval_same = .true.
640  self%index_shift = (self%spline_degree-2)/2
641  else
642  self%eval_same = .false.
643  self%index_shift = (self%spline_degree-1)/2
644  end if
645  else
646  if ( modulo( self%spline_degree, 2 ) == 0 ) then
647  self%eval_same = .false.
648  self%index_shift = (self%spline_degree)/2
649  else
650  self%eval_same = .true.
651  self%index_shift = (self%spline_degree-1)/2
652  end if
653  end if
654  end if
655 
656  ! normalized Gauss Legendre points and weights
657  self%n_quad_points = (self%spline_degree+2)/2
658  self%quad_xw = sll_f_gauss_legendre_points_and_weights(self%n_quad_points)
659 !!$
660 !!$ ! For nonlinear disgradE
661 !!$ self%n_quadp1_points = (self%spline_degree+3)/2
662 !!$ allocate( self%quadp1_xw(2, self%n_quadp1_points) )
663 !!$ allocate( self%spline_pol_val(self%n_span) )
664 !!$ ! normalized Gauss Legendre points and weights
665 !!$ self%quadp1_xw = sll_f_gauss_legendre_points_and_weights(self%n_quadp1_points)
666 !!$
667 !!$
668 !!$ ! For smoothed evaluate and add_charge
669 !!$ allocate( self%quads1_xw(2, self%n_quadp1_points) )
670 !!$ ! normalized Gauss Legendre points and weights to [0,1]
671 !!$ self%quads1_xw = self%quadp1_xw
672 !!$ self%quads1_xw(1,:) = (self%quads1_xw(1,:) + 1.0_f64)*0.5_f64
673 !!$ self%quads1_xw(2,:) = self%quads1_xw(2,:)*0.5_f64
674 !!$
675 !!$
676 !!$
677 !!$ ! For smoothed add_current
678 !!$ self%n_quads2_points = (self%spline_degree+4)/2
679 !!$ allocate( self%quads2_xw(2, self%n_quads2_points) )
680 !!$ ! normalized Gauss Legendre points and weights to [0,1]
681 !!$ self%quads2_xw = sll_f_gauss_legendre_points_and_weights(self%n_quads2_points)
682 !!$ self%quads2_xw(1,:) = (self%quads2_xw(1,:) + 1.0_f64)*0.5_f64
683 !!$ self%quads2_xw(2,:) = self%quads2_xw(2,:)*0.5_f64
684 
685  end subroutine init_spline_strong_1d
686 
688  subroutine free_spline_strong_1d(self)
689  class(sll_t_particle_mesh_coupling_spline_strong_1d), intent( inout ) :: self
690 
691  deallocate(self%spline_val)
692  deallocate(self%spline_val_more)
693  deallocate(self%quad_xw)
694  deallocate(self%spline_val_prim)
695  deallocate(self%spline_val_prim_old)
696  deallocate(self%spline_val_prim_new)
697  if (self%spline_degree < 6 ) then
698  call sll_s_spline_pp_free_1d(self%spline_pp)
699  end if
700 
701  end subroutine free_spline_strong_1d
702 
703 !!!!!!!!!!!! FROM HERE NOT IMPLEMENTED !!!!!!!!!!!!!
704 
705 
706  !---------------------------------------------------------------------------!
708  subroutine add_particle_mass_spline_strong_1d(self, position, marker_charge, particle_mass)
709  class( sll_t_particle_mesh_coupling_spline_strong_1d ), intent(inout) :: self
710  sll_real64, intent( in ) :: position(self%dim)
711  sll_real64, intent( in ) :: marker_charge
712  sll_real64, intent( inout ) :: particle_mass(:,:)
713 
714 
715  sll_error("add_particle_mass", "Not implemented.")
716 
717  end subroutine add_particle_mass_spline_strong_1d
718 
719  !---------------------------------------------------------------------------!
721  subroutine evaluate_multiple_spline_strong_1d(self, position, components, field_dofs, field_value)
722  class(sll_t_particle_mesh_coupling_spline_strong_1d), intent( inout ) :: self
723  sll_real64, intent( in ) :: position(self%dim)
724  sll_int32, intent(in) :: components(:)
725  sll_real64, intent( in ) :: field_dofs(:,:)
726  sll_real64, intent(out) :: field_value(:)
727 
728  sll_error("evaluate_multiple", "Not implemented.")
729 
730  end subroutine evaluate_multiple_spline_strong_1d
731 
732 
733  subroutine evaluate_field_single_spline_strong_pp_1d(self, position, field_dofs_pp, field_value)
734  class(sll_t_particle_mesh_coupling_spline_strong_1d), intent( inout ) :: self
735  sll_real64, intent( in ) :: position(self%dim)
736  sll_real64, intent(in) :: field_dofs_pp(:,:)
737  sll_real64, intent( out ) :: field_value
738 
739  sll_error("evaluate_pp", "Not implemented.")
740 
741  end subroutine evaluate_field_single_spline_strong_pp_1d
742 
743 
745  subroutine add_current_spline_strong_1d( self, position_old, position_new, marker_charge, j_dofs )
746  class( sll_t_particle_mesh_coupling_spline_strong_1d ), intent(inout) :: self
747  sll_real64, intent( in ) :: position_old(self%dim)
748  sll_real64, intent( in ) :: position_new(self%dim)
749  sll_real64, intent( in ) :: marker_charge
750  sll_real64, intent( inout ) :: j_dofs(self%n_dofs)
751 
752 
753  sll_error("add_current", "Not implemented.")
754 
755  end subroutine add_current_spline_strong_1d
756 
758  subroutine add_current_update_v_pp_spline_strong_1d (self, position_old, position_new, marker_charge, qoverm, bfield_dofs, vi, j_dofs)
759  class(sll_t_particle_mesh_coupling_spline_strong_1d), intent(inout) :: self
760  sll_real64, intent(in) :: position_old(self%dim)
761  sll_real64, intent(in) :: position_new(self%dim)
762  sll_real64, intent(in) :: marker_charge
763  sll_real64, intent(in) :: qoverm
764  sll_real64, intent(in) :: bfield_dofs(self%n_dofs)
765  sll_real64, intent(inout) :: vi(:)
766  sll_real64, intent(inout) :: j_dofs(self%n_dofs)
767 
768 
769  sll_error("add_current_update_v_pp", "Not implemented.")
770 
771  end subroutine add_current_update_v_pp_spline_strong_1d
772 
773  !----------------------------------------------------------------------------------
774 
775  subroutine sll_s_new_particle_mesh_coupling_spline_strong_1d( smoother, domain, n_cells, spline_degree, integ, eval_grid_points )
776  class( sll_c_particle_mesh_coupling_1d), allocatable, intent(out):: smoother
777  sll_int32, intent(in) :: n_cells
778  sll_real64, intent(in) :: domain(2)
779  sll_int32, intent(in) :: spline_degree
780  logical, intent( in ) :: integ
781  logical, intent( in ) :: eval_grid_points
782 
783  !local variables
784  sll_int32 :: ierr
785 
786 
787  allocate( sll_t_particle_mesh_coupling_spline_strong_1d :: smoother , stat=ierr)
788  sll_assert( ierr == 0)
789 
790  select type( smoother )
791  type is ( sll_t_particle_mesh_coupling_spline_strong_1d )
792  call smoother%init( domain, n_cells, spline_degree, integ, eval_grid_points )
793  end select
794 
795  end subroutine sll_s_new_particle_mesh_coupling_spline_strong_1d
796 
797 end module sll_m_particle_mesh_coupling_spline_strong_1d
sll_m_particle_mesh_coupling_base_1d::add_current_evaluate
Definition: sll_m_particle_mesh_coupling_base_1d.F90:115
sll_m_particle_mesh_coupling_base_1d::add_current
Definition: sll_m_particle_mesh_coupling_base_1d.F90:101
sll_m_gauss_legendre_integration
Gauss-Legendre integration.
Definition: sll_m_gauss_legendre_integration.F90:8
sll_m_gauss_legendre_integration::sll_f_gauss_legendre_points_and_weights
real(kind=f64) function, dimension(2, 1:npoints), public sll_f_gauss_legendre_points_and_weights(npoints, a, b)
Returns a 2d array of size (2,npoints) containing gauss-legendre points and weights in the interval [...
Definition: sll_m_gauss_legendre_integration.F90:265
sll_m_low_level_bsplines
Low level arbitrary degree splines.
Definition: sll_m_low_level_bsplines.F90:27
sll_m_particle_group_base
Definition: sll_m_particle_group_base.F90:1
sll_m_particle_mesh_coupling_base_1d
Base class for kernel smoothers for accumulation and field evaluation in PIC.
Definition: sll_m_particle_mesh_coupling_base_1d.F90:5
sll_m_particle_mesh_coupling_base_1d::sll_p_collocation
integer(kind=i32), parameter, public sll_p_collocation
Definition: sll_m_particle_mesh_coupling_base_1d.F90:22
sll_m_particle_mesh_coupling_base_1d::sll_p_galerkin
integer(kind=i32), parameter, public sll_p_galerkin
Definition: sll_m_particle_mesh_coupling_base_1d.F90:21
sll_m_particle_mesh_coupling_spline_strong_1d
Kernel smoother for 1d with splines of arbitrary degree placed on a uniform mesh. This version is for...
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:6
sll_m_particle_mesh_coupling_spline_strong_1d::add_particle_mass_spline_strong_1d
subroutine add_particle_mass_spline_strong_1d(self, position, marker_charge, particle_mass)
Add charge of one particle.
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:709
sll_m_particle_mesh_coupling_spline_strong_1d::evaluate_multiple_spline_strong_1d
subroutine evaluate_multiple_spline_strong_1d(self, position, components, field_dofs, field_value)
Evaluate several fields at position position.
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:722
sll_m_particle_mesh_coupling_spline_strong_1d::add_current_evaluate_spline_strong_1d
subroutine add_current_evaluate_spline_strong_1d(self, position_old, position_new, marker_charge, vbar, field_dofs, j_dofs, field)
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:481
sll_m_particle_mesh_coupling_spline_strong_1d::add_charge_single_spline_strong_1d
subroutine add_charge_single_spline_strong_1d(self, position, marker_charge, rho_dofs)
Add charge of one particle.
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:129
sll_m_particle_mesh_coupling_spline_strong_1d::evaluate_field_single_spline_strong_pp_1d
subroutine evaluate_field_single_spline_strong_pp_1d(self, position, field_dofs_pp, field_value)
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:734
sll_m_particle_mesh_coupling_spline_strong_1d::add_current_update_v_spline_strong_1d
subroutine add_current_update_v_spline_strong_1d(self, position_old, position_new, marker_charge, qoverm, bfield_dofs, vi, j_dofs)
Add current for one particle and update v (according to H_p1 part in Hamiltonian splitting)
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:256
sll_m_particle_mesh_coupling_spline_strong_1d::add_current_update_v_pp_spline_strong_1d
subroutine add_current_update_v_pp_spline_strong_1d(self, position_old, position_new, marker_charge, qoverm, bfield_dofs, vi, j_dofs)
Add current for one particle and update v (according to H_p1 part in Hamiltonian splitting)
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:759
sll_m_particle_mesh_coupling_spline_strong_1d::add_current_update_v_spline_strong_1d_quadrature
subroutine add_current_update_v_spline_strong_1d_quadrature(self, position_old, position_new, marker_charge, qoverm, bfield_dofs, vi, j_dofs)
Add current for one particle and update v (according to H_p1 part in Hamiltonian splitting)....
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:172
sll_m_particle_mesh_coupling_spline_strong_1d::add_current_evaluate_spline_strong_1d_quadrature
subroutine add_current_evaluate_spline_strong_1d_quadrature(self, position_old, position_new, marker_charge, vbar, field_dofs, j_dofs, field)
Add current for one particle and update v (according to H_p1 part in Hamiltonian splitting)....
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:393
sll_m_particle_mesh_coupling_spline_strong_1d::helper_normalized_position
subroutine helper_normalized_position(self, position, index, axi)
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:83
sll_m_particle_mesh_coupling_spline_strong_1d::add_current_spline_strong_1d
subroutine add_current_spline_strong_1d(self, position_old, position_new, marker_charge, j_dofs)
Add current with integration over x.
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:746
sll_m_particle_mesh_coupling_spline_strong_1d::free_spline_strong_1d
subroutine free_spline_strong_1d(self)
Destructor.
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:689
sll_m_particle_mesh_coupling_spline_strong_1d::current_eval_local
subroutine current_eval_local(self, upper, lower, box, marker_charge, field, j_dofs, field_dofs)
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:435
sll_m_particle_mesh_coupling_spline_strong_1d::current_v_local
subroutine current_v_local(self, upper, lower, box, marker_charge, qoverm, vi, j_dofs, bfield_dofs)
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:210
sll_m_particle_mesh_coupling_spline_strong_1d::sll_s_new_particle_mesh_coupling_spline_strong_1d
subroutine, public sll_s_new_particle_mesh_coupling_spline_strong_1d(smoother, domain, n_cells, spline_degree, integ, eval_grid_points)
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:776
sll_m_particle_mesh_coupling_spline_strong_1d::init_spline_strong_1d
subroutine init_spline_strong_1d(self, domain, n_cells, spline_degree, integ, eval_grid_points)
Constructor.
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:563
sll_m_particle_mesh_coupling_spline_strong_1d::evaluate_field_single_spline_strong_1d
subroutine evaluate_field_single_spline_strong_1d(self, position, field_dofs, field_value)
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:344
sll_m_splines_pp
Splines in pp form.
Definition: sll_m_splines_pp.F90:8
sll_m_splines_pp::sll_s_spline_pp_init_1d
subroutine, public sll_s_spline_pp_init_1d(self, degree, n_cells, boundary)
Initialize sll_t_spline_pp_1d object (Set poly_coeffs depending on degree)
Definition: sll_m_splines_pp.F90:1387
sll_m_splines_pp::sll_s_spline_pp_free_1d
subroutine, public sll_s_spline_pp_free_1d(self)
Destructor 1d.
Definition: sll_m_splines_pp.F90:1685
sll_m_splines_pp::sll_f_spline_pp_horner_1d
real(kind=f64) function, public sll_f_spline_pp_horner_1d(degree, pp_coeffs, x, index)
Perform a 1d hornerschema on the pp_coeffs at index.
Definition: sll_m_splines_pp.F90:1175
sll_m_particle_group_base::sll_c_particle_group_base
Definition: sll_m_particle_group_base.F90:42
sll_m_particle_mesh_coupling_base_1d::sll_c_particle_mesh_coupling_1d
Basic type of a kernel smoother used for PIC simulations.
Definition: sll_m_particle_mesh_coupling_base_1d.F90:25
sll_m_particle_mesh_coupling_spline_strong_1d::sll_t_particle_mesh_coupling_spline_strong_1d
Basic type of a kernel smoother used for PIC simulations.
Definition: sll_m_particle_mesh_coupling_spline_strong_1d.F90:38
sll_m_splines_pp::sll_t_spline_pp_1d
arbitrary degree 1d spline
Definition: sll_m_splines_pp.F90:93
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:40 for Semi-Lagrangian Library by   doxygen 1.9.1