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_1d.F90
Go to the documentation of this file.
1 
6 
7  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8 #include "sll_assert.h"
9 #include "sll_memory.h"
10 #include "sll_working_precision.h"
11 #include "sll_errors.h"
12 
15 
16  use sll_m_low_level_bsplines, only: &
17  sll_s_uniform_bsplines_eval_basis, &
18  sll_s_uniform_bsplines_eval_deriv
19 
24 
25  use sll_m_particle_group_base, only: &
27 
28  use sll_m_splines_pp, only: &
34 
35  implicit none
36 
37  public :: &
41 
42  private
43  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
44 
47  type(sll_t_spline_pp_1d) :: spline_pp
48 
49  ! Information about the particles
50  sll_int32 :: no_particles
51  sll_int32 :: n_span
52  sll_real64 :: scaling
53  sll_int32 :: n_quad_points
54  sll_int32 :: n_quadp1_points
55  sll_int32 :: n_quadp2_points
56  sll_int32 :: n_quads2_points
57 
58  sll_real64, allocatable :: spline_val(:)
59  sll_real64, allocatable :: spline_pol_val(:)
60  sll_real64, allocatable :: spline_val_more(:)
61  sll_real64, allocatable :: quad_xw(:,:)
62  sll_real64, allocatable :: quadp1_xw(:,:)
63  sll_real64, allocatable :: quadp2_xw(:,:)
64  sll_real64, allocatable :: quads1_xw(:,:)
65  sll_real64, allocatable :: quads2_xw(:,:)
66 
67  sll_int32 :: n_quad_points_line
68  sll_real64, allocatable :: quad_xw_line(:,:)
69 
70  contains
71  procedure :: add_charge => add_charge_single_spline_1d
72  procedure :: add_particle_mass => add_particle_mass_spline_1d
73  procedure :: add_particle_mass_full => add_particle_mass_spline_1d_full
74  procedure :: evaluate => evaluate_field_single_spline_1d
75  procedure :: evaluate_multiple => evaluate_multiple_spline_1d
76  procedure :: add_current_update_v => add_current_update_v_spline_pp_1d
78 
79  procedure :: init => init_spline_1d
80  procedure :: free => free_spline_1d
81 
83  procedure :: evaluate_int => evaluate_int_spline_1d
84  procedure :: evaluate_int_quad => evaluate_int_quad_spline_1d
85  procedure :: evaluate_int_linear_quad => evaluate_int_linear_quad_spline_1d
86 
87  procedure :: evaluate_simple => evaluate_field_single_spline_1d
88 
89  procedure :: evaluate_int_subnewton => evaluate_int_spline_1d_subnewton
90  procedure :: evaluate_int_linear_quad_subnewton => evaluate_int_linear_quad_subnewton_spline_1d
91  procedure :: add_current_split => add_current_split_spline_1d
92  procedure :: add_charge_int => add_charge_int_spline_1d
93 
94  procedure :: update_jv
95 
97 
98 contains
99 
100 
101  !---------------------------------------------------------------------------!
103  subroutine add_charge_single_spline_1d(self, position, marker_charge, rho_dofs)
104  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
105  sll_real64, intent( in ) :: position(self%dim)
106  sll_real64, intent( in ) :: marker_charge
107  sll_real64, intent( inout ) :: rho_dofs(self%n_cells)
108  !local variables
109  sll_int32 :: i1
110  sll_int32 :: index1d, box
111  sll_real64 :: xi
112 
113  call convert_x_to_xbox( self, position, xi, box )
114  box = box - self%spline_degree
115  call sll_s_uniform_bsplines_eval_basis(self%spline_degree, xi, self%spline_val)
116  !alternative version with horner and pp-form
117  !call sll_s_spline_pp_horner_m_1d(self%spline_pp, self%spline_val, self%spline_degree, xi)
118 
119  do i1 = 1, self%n_span
120  index1d = modulo(box+i1-2,self%n_cells)+1
121  rho_dofs(index1d) = rho_dofs(index1d) +&
122  (marker_charge * self%spline_val(i1)* self%scaling)
123  end do
124 
125  end subroutine add_charge_single_spline_1d
126 
127 
128  !---------------------------------------------------------------------------!
130  subroutine add_particle_mass_spline_1d(self, position, marker_charge, particle_mass)
131  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
132  sll_real64, intent( in ) :: position(self%dim)
133  sll_real64, intent( in ) :: marker_charge
134  sll_real64, intent( inout ) :: particle_mass(:,:)
135 
136  !local variables
137  sll_int32 :: i1, column
138  sll_int32 :: index1d, box
139  sll_real64 :: xi
140 
141  sll_assert( size(particle_mass,1) == self%n_span )
142  sll_assert( size(particle_mass,2) == self%n_cells )
143 
144  call convert_x_to_xbox( self, position, xi, box )
145  box = box - self%spline_degree
146  call sll_s_uniform_bsplines_eval_basis(self%spline_degree, xi, self%spline_val)
147  !alternative version with horner and pp-form
148  !call sll_s_spline_pp_horner_m_1d(self%spline_pp, self%spline_val, self%spline_degree, xi)
149 
150  do i1 = 1, self%n_span
151  index1d = modulo(box+i1-2,self%n_cells)+1
152  do column = i1, self%n_span
153  particle_mass(column-i1+1, index1d) = particle_mass( column-i1+1, index1d) + &
154  marker_charge * self%spline_val(i1) * self%spline_val(column)
155  ! Note: No scaling since Galerkin function
156  end do
157  end do
158 
159  end subroutine add_particle_mass_spline_1d
160 
161  !----------------------------------------------------------------------!
163  subroutine add_particle_mass_spline_1d_full(self, position, marker_charge, particle_mass)
164  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
165  sll_real64, intent( in ) :: position(self%dim)
166  sll_real64, intent( in ) :: marker_charge
167  sll_real64, intent( inout ) :: particle_mass(:,:)
168 
169  !local variables
170  sll_int32 :: i1, column, ind
171  sll_int32 :: index1d, box
172  sll_real64 :: xi
173 
174  sll_assert( size(particle_mass,1) == 2*self%spline_degree+1 )
175  sll_assert( size(particle_mass,2) == self%n_cells )
176 
177  call convert_x_to_xbox( self, position, xi, box )
178  box = box - self%spline_degree
179  call sll_s_uniform_bsplines_eval_basis(self%spline_degree, xi, self%spline_val)
180  !alternative version with horner and pp-form
181  !call sll_s_spline_pp_horner_m_1d(self%spline_pp, self%spline_val, self%spline_degree, xi)
182  do i1 = 1, self%n_span
183  index1d = modulo(box+i1-2,self%n_cells)+1
184  ind=1+(self%n_span-i1)
185  do column = 1, self%n_span
186  particle_mass(ind, index1d) = particle_mass( ind, index1d) + &
187  marker_charge * self%spline_val(i1) * self%spline_val(column)
188  ind = ind+1
189  end do
190  end do
191 
192  end subroutine add_particle_mass_spline_1d_full
193 
194 
195  !---------------------------------------------------------------------------!
197  subroutine evaluate_field_single_spline_1d(self, position, field_dofs, field_value)
198  class(sll_t_particle_mesh_coupling_spline_1d), intent( inout ) :: self
199  sll_real64, intent( in ) :: position(self%dim)
200  sll_real64, intent( in ) :: field_dofs(self%n_cells)
201  sll_real64, intent( out ) :: field_value
202 
203  !local variables
204  sll_int32 :: i1
205  sll_int32 :: index1d, box
206  sll_real64 :: xi
207 
208  call convert_x_to_xbox( self, position, xi, box )
209  box = box - self%spline_degree
210  ! Version based on the uniform bspline functions
211  call sll_s_uniform_bsplines_eval_basis(self%spline_degree, xi, self%spline_val)
212  !alternative version with horner and pp-form
213  !call sll_s_spline_pp_horner_m_1d(self%spline_pp, self%spline_val, self%spline_degree, xi)
214 
215  field_value = 0.0_f64
216  do i1 = 1, self%n_span
217  index1d = modulo(box+i1-2, self%n_cells)+1
218  field_value = field_value + &
219  field_dofs(index1d) * &
220  self%spline_val(i1)
221  end do
222 
223  end subroutine evaluate_field_single_spline_1d
224 
225 
226  !---------------------------------------------------------------------------!
228  subroutine evaluate_multiple_spline_1d(self, position, components, field_dofs, field_value)
229  class(sll_t_particle_mesh_coupling_spline_1d), intent( inout ) :: self
230  sll_real64, intent( in ) :: position(self%dim)
231  sll_int32, intent(in) :: components(:)
232  sll_real64, intent( in ) :: field_dofs(:,:)
233  sll_real64, intent(out) :: field_value(:)
234 
235  !local variables
236  sll_int32 :: i1
237  sll_int32 :: index1d, box
238  sll_real64 :: xi
239 
240  sll_assert( size(field_dofs,1) == self%n_cells )
241  sll_assert( size(field_dofs,2) == size(field_value) )
242 
243  call convert_x_to_xbox( self, position, xi, box )
244 
245  box = box - self%spline_degree
246  call sll_s_uniform_bsplines_eval_basis(self%spline_degree, xi, self%spline_val)
247  !alternative version with horner and pp-form
248  !call sll_s_spline_pp_horner_m_1d(self%spline_pp, self%spline_val, self%spline_degree, xi)
249 
250  field_value = 0.0_f64
251  do i1 = 1, self%n_span
252  index1d = modulo(box+i1-2, self%n_cells)+1
253  field_value = field_value + &
254  field_dofs(index1d,components) * &
255  self%spline_val(i1)
256  end do
257 
258  end subroutine evaluate_multiple_spline_1d
259 
260 
262  subroutine add_current_spline_1d( self, position_old, position_new, marker_charge, j_dofs )
263  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
264  sll_real64, intent( in ) :: position_old(self%dim)
265  sll_real64, intent( in ) :: position_new(self%dim)
266  sll_real64, intent( in ) :: marker_charge
267  sll_real64, intent( inout ) :: j_dofs(self%n_cells)
268 
269  sll_real64 :: dx_new, dx_old
270  sll_int32 :: box_new, box_old, i, index, j, ind
271 
272  call convert_x_to_xbox( self, position_new, dx_new, box_new )
273  call convert_x_to_xbox( self, position_old, dx_old, box_old )
274 
275  !print*, dx_old, box_old, dx_new, box_new
276 
277  do i=1, self%n_span
278  self%spline_val(i) = sll_f_spline_pp_horner_1d(self%spline_degree+1, self%spline_pp%poly_coeffs_fpa, dx_old, i) &
279  * self%delta_x
280  self%spline_val_more(i) = sll_f_spline_pp_horner_1d(self%spline_degree+1, self%spline_pp%poly_coeffs_fpa, dx_new, i) &
281  * self%delta_x
282  end do
283 
284  ! New version to avoid the cancellation
285  if (box_old<box_new) then
286  i=0
287  do ind = box_old-self%spline_degree, min(box_new-self%spline_degree-1, box_old )
288  index = modulo(ind-1,self%n_cells)+1
289  i = i+1
290  j_dofs(index) = j_dofs(index) + (marker_charge * &
291  (self%delta_x-self%spline_val(i)))
292  end do
293  j=0
294  do ind = box_new-self%spline_degree, box_old
295  j = j+1
296  i = i+1
297  index = modulo(ind-1,self%n_cells)+1
298  j_dofs(index) = j_dofs(index) + marker_charge * &
299  ((self%spline_val_more(j) - self%spline_val(i)))
300  end do
301  do ind = max(box_old+1, box_new-self%spline_degree), box_new
302  j = j+1
303  index = modulo(ind-1,self%n_cells)+1
304  j_dofs(index) = j_dofs(index) + (marker_charge * &
305  (self%spline_val_more(j) ))
306  end do
307 
308  do ind = box_old+1, box_new-self%spline_degree-1
309  index = modulo(ind-1, self%n_cells)+1
310  j_dofs(index) = j_dofs(index) + marker_charge * self%delta_x
311  end do
312  else
313 
314  j=0
315  do ind = box_new-self%spline_degree, min(box_old-self%spline_degree-1, box_new)
316  index = modulo(ind-1,self%n_cells)+1
317  j = j+1
318  j_dofs(index) = j_dofs(index) + (marker_charge * &
319  (-self%delta_x+self%spline_val_more(j)))
320  end do
321  i=0
322  do ind = box_old-self%spline_degree, box_new
323  j = j+1
324  i = i+1
325  index = modulo(ind-1,self%n_cells)+1
326  j_dofs(index) = j_dofs(index) + (marker_charge * &
327  (self%spline_val_more(j) - self%spline_val(i)))
328  end do
329  do ind = max(box_new+1, box_old-self%spline_degree), box_old
330  i = i+1
331  index = modulo(ind-1,self%n_cells)+1
332  j_dofs(index) = j_dofs(index) + (marker_charge * &
333  (-self%spline_val(i) ))
334  end do
335 
336  do ind = box_new+1, box_old-self%spline_degree-1
337  index = modulo(ind-1, self%n_cells)+1
338  j_dofs(index) = j_dofs(index) - marker_charge * self%delta_x
339  end do
340  end if
341 
342  end subroutine add_current_spline_1d
343 
344 
346  subroutine add_current_evaluate_spline_1d(self, position_old, position_new, marker_charge, vbar, field_dofs, j_dofs, field)
347  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
348  sll_real64, intent( in ) :: position_old(self%dim)
349  sll_real64, intent( in ) :: position_new(self%dim)
350  sll_real64, intent( in ) :: marker_charge
351  sll_real64, intent( in ) :: vbar
352  sll_real64, intent( in ) :: field_dofs(self%n_dofs)
353  sll_real64, intent( inout ) :: j_dofs(self%n_dofs)
354  sll_real64, intent( out ) :: field
355 
356  sll_real64 :: dx_new, dx_old
357  sll_int32 :: box_new, box_old, i, index, ind, j
358 
359  call convert_x_to_xbox( self, position_new, dx_new, box_new )
360  call convert_x_to_xbox( self, position_old, dx_old, box_old )
361 
362 
363  do i=1, self%n_span
364  self%spline_val(i) = sll_f_spline_pp_horner_1d(self%spline_degree+1, self%spline_pp%poly_coeffs_fpa, dx_old, i) &
365  * self%delta_x
366  self%spline_val_more(i) = sll_f_spline_pp_horner_1d(self%spline_degree+1, self%spline_pp%poly_coeffs_fpa, dx_new, i) &
367  * self%delta_x
368  end do
369 
370  field = 0.0_f64
371 
372 
373  ! New version to avoid the cancellation
374  if (box_old<box_new) then
375  i=0
376  do ind = box_old-self%spline_degree, min(box_new-self%spline_degree-1, box_old)
377  index = modulo(ind-1,self%n_cells)+1
378  i = i+1
379  j_dofs(index) = j_dofs(index) + (marker_charge * &
380  (self%delta_x-self%spline_val(i)))
381  field = field + ((self%delta_x - self%spline_val(i))) * field_dofs(index)
382  end do
383  j=0
384  do ind = box_new-self%spline_degree, box_old
385  j = j+1
386  i = i+1
387  index = modulo(ind-1,self%n_cells)+1
388  j_dofs(index) = j_dofs(index) + marker_charge * &
389  ((self%spline_val_more(j) - self%spline_val(i)))
390  field = field + ((self%spline_val_more(j) - self%spline_val(i))) * field_dofs(index)
391  end do
392 
393  do ind = max(box_old+1, box_new-self%spline_degree), box_new
394  j = j+1
395  index = modulo(ind-1,self%n_cells)+1
396  j_dofs(index) = j_dofs(index) + (marker_charge * &
397  (self%spline_val_more(j) ))
398  field = field + (self%spline_val_more(j) ) * field_dofs(index)
399  end do
400 
401  do ind = box_old+1, box_new-self%spline_degree-1
402  index = modulo(ind-1, self%n_cells)+1
403  j_dofs(index) = j_dofs(index) + marker_charge * self%delta_x
404  field = field + self%delta_x * field_dofs(index)
405  end do
406  else
407 
408  j=0
409  do ind = box_new-self%spline_degree, min(box_old-self%spline_degree-1, box_new)
410  index = modulo(ind-1,self%n_cells)+1
411  j = j+1
412  j_dofs(index) = j_dofs(index) + (marker_charge * &
413  (-self%delta_x+self%spline_val_more(j)))
414  field = field + (-self%delta_x+self%spline_val_more(j)) * field_dofs(index)
415  end do
416  i=0
417  do ind = box_old-self%spline_degree, box_new
418  j = j+1
419  i = i+1
420  index = modulo(ind-1,self%n_cells)+1
421  j_dofs(index) = j_dofs(index) + (marker_charge * &
422  (self%spline_val_more(j) - self%spline_val(i)))
423  field = field + (self%spline_val_more(j) - self%spline_val(i)) * field_dofs(index)
424  end do
425  do ind = max(box_new+1, box_old-self%spline_degree), box_old
426  i = i+1
427  index = modulo(ind-1,self%n_cells)+1
428  j_dofs(index) = j_dofs(index) + (marker_charge * &
429  (-self%spline_val(i) ))
430  field = field + (-self%spline_val(i) ) * field_dofs(index)
431  end do
432 
433  do ind = box_new+1, box_old-self%spline_degree-1
434  index = modulo(ind-1, self%n_cells)+1
435  j_dofs(index) = j_dofs(index) - marker_charge * self%delta_x
436  field = field - self%delta_x * field_dofs(index)
437  end do
438 
439  end if
440 
441  field = field/vbar
442 
443  end subroutine add_current_evaluate_spline_1d
444 
446  subroutine evaluate_int_spline_1d(self, position_old, position_new, vbar, field_dofs, field)
447  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
448  sll_real64, intent( in ) :: position_old(self%dim)
449  sll_real64, intent( in ) :: position_new(self%dim)
450  sll_real64, intent( in ) :: vbar
451  sll_real64, intent( in ) :: field_dofs(self%n_cells)
452  sll_real64, intent( out ) :: field
453 
454  sll_real64 :: dx_new, dx_old
455  sll_int32 :: box_new, box_old, i, index, ind, j
456 
457  call convert_x_to_xbox( self, position_new, dx_new, box_new )
458  call convert_x_to_xbox( self, position_old, dx_old, box_old )
459 
460 
461  do i=1, self%n_span
462  self%spline_val(i) = sll_f_spline_pp_horner_1d(self%spline_degree+1, self%spline_pp%poly_coeffs_fpa, dx_old, i) &
463  * self%delta_x
464  self%spline_val_more(i) = sll_f_spline_pp_horner_1d(self%spline_degree+1, self%spline_pp%poly_coeffs_fpa, dx_new, i) &
465  * self%delta_x
466  end do
467 
468  field = 0.0_f64
469 
470 
471  ! New version to avoid the cancellation
472  if (box_old<box_new) then
473  i=0
474  do ind = box_old-self%spline_degree, min(box_new-self%spline_degree-1, box_old)
475  index = modulo(ind-1,self%n_cells)+1
476  i = i+1
477  field = field + ((self%delta_x - self%spline_val(i))) * field_dofs(index)
478  end do
479  j=0
480  do ind = box_new-self%spline_degree, box_old
481  j = j+1
482  i = i+1
483  index = modulo(ind-1,self%n_cells)+1
484  field = field + ((self%spline_val_more(j) - self%spline_val(i))) * field_dofs(index)
485  end do
486  do ind = max(box_old+1, box_new-self%spline_degree), box_new
487  j = j+1
488  index = modulo(ind-1,self%n_cells)+1
489  field = field + (self%spline_val_more(j) ) * field_dofs(index)
490  end do
491  do ind = box_old+1, box_new-self%spline_degree-1
492  index = modulo(ind-1, self%n_cells)+1
493  field = field + self%delta_x * field_dofs(index)
494  end do
495  else
496 
497  j=0
498  do ind = box_new-self%spline_degree, min(box_old-self%spline_degree-1, box_new)
499  index = modulo(ind-1,self%n_cells)+1
500  j = j+1
501  field = field + (-self%delta_x+self%spline_val_more(j)) * field_dofs(index)
502  end do
503  i=0
504  do ind = box_old-self%spline_degree, box_new
505  j = j+1
506  i = i+1
507  index = modulo(ind-1,self%n_cells)+1
508  field = field + (self%spline_val_more(j) - self%spline_val(i)) * field_dofs(index)
509  end do
510  do ind = max(box_new+1,box_old-self%spline_degree), box_old
511  i = i+1
512  index = modulo(ind-1,self%n_cells)+1
513  field = field + (-self%spline_val(i) ) * field_dofs(index)
514  end do
515 
516  do ind = box_new+1, box_old-self%spline_degree-1
517  index = modulo(ind-1, self%n_cells)+1
518  field = field - self%delta_x * field_dofs(index)
519  end do
520 
521  end if
522 
523  field = field/vbar
524 
525  end subroutine evaluate_int_spline_1d
526 
528  subroutine evaluate_int_quad_spline_1d(self, position_old, position_new, field_dofs, field)
529  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
530  sll_real64, intent( in ) :: position_old(self%dim)
531  sll_real64, intent( in ) :: position_new(self%dim)
532  sll_real64, intent( in ) :: field_dofs(self%n_cells)
533  sll_real64, intent( out ) :: field
534  !local variables
535  sll_real64 :: dx_new, dx_old
536  sll_int32 :: box_new, box_old, ind
537  sll_real64 :: rinterval, lowert, uppert
538 
539  call convert_x_to_xbox( self, position_new, dx_new, box_new )
540  call convert_x_to_xbox( self, position_old, dx_old, box_old )
541 
542  rinterval = self%delta_x/(position_new(1)-position_old(1))
543 
544  field = 0.0_f64
545 
546  if ( box_old == box_new ) then
547  call update_v_partial(self, dx_new, dx_old, 1.0_f64, 0.0_f64, box_old, field_dofs, field )
548  elseif ( box_old < box_new ) then
549  lowert = 0.0_f64
550  uppert = (1.0_f64-dx_old)*rinterval
551  call update_v_partial( self, 1.0_f64, dx_old, uppert, lowert, box_old, field_dofs, field )
552  do ind = box_old+1, box_new-1
553  lowert = uppert
554  uppert = uppert + rinterval
555  call update_v_partial( self, 1.0_f64, 0.0_f64, uppert, lowert, ind, field_dofs, field )
556  end do
557  call update_v_partial( self, dx_new, 0.0_f64, 1.0_f64, uppert, box_new, field_dofs, field )
558  else
559  lowert = 1.0_f64
560  uppert = 1.0_f64 + (rinterval-dx_new*rinterval)
561  call update_v_partial( self, 1.0_f64, dx_new, uppert, lowert, box_new, field_dofs, field )
562  do ind = box_new+1, box_old-1
563  lowert = uppert
564  uppert = lowert + rinterval
565  call update_v_partial( self, 1.0_f64, 0.0_f64, uppert, lowert, ind, field_dofs, field )
566  end do
567  call update_v_partial( self, dx_old, 0.0_f64, 0.0_f64, uppert, box_old, field_dofs, field )
568  field = - field ! Account for the fact that we have integrated the wrong direction
569  end if
570 
571  end subroutine evaluate_int_quad_spline_1d
572 
573  ! Evaluation of the integral needed for Newton's method in subcycling
574  subroutine evaluate_int_spline_1d_subnewton(self, position_old, position_new, field_dofs, field)
575  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
576  sll_real64, intent( in ) :: position_old(self%dim)
577  sll_real64, intent( in ) :: position_new(self%dim)
578  sll_real64, intent( in ) :: field_dofs(self%n_cells)
579  sll_real64, intent( out ) :: field
580  !local variables
581  sll_real64 :: dx_new, dx_old
582  sll_int32 :: box_new, box_old, ind
583  sll_real64 :: rinterval, lowert, uppert
584 
585  call convert_x_to_xbox( self, position_new, dx_new, box_new )
586  call convert_x_to_xbox( self, position_old, dx_old, box_old )
587 
588  rinterval = self%delta_x/(position_new(1)-position_old(1))
589 
590  field = 0.0_f64
591 
592  if ( box_old == box_new ) then
593  call update_v_partial_newton(self, dx_new, dx_old, 1.0_f64, 0.0_f64, box_old, field_dofs, field )
594  elseif ( box_old < box_new ) then
595  lowert = 0.0_f64
596  uppert = (1.0_f64-dx_old)*rinterval
597  call update_v_partial_newton( self, 1.0_f64, dx_old, uppert, lowert, box_old, field_dofs, field )
598  do ind = box_old+1, box_new-1
599  lowert = uppert
600  uppert = uppert + rinterval
601  call update_v_partial_newton( self, 1.0_f64, 0.0_f64, uppert, lowert, ind, field_dofs, field )
602  end do
603  call update_v_partial_newton( self, dx_new, 0.0_f64, 1.0_f64, uppert, box_new, field_dofs, field )
604  else
605  lowert = 1.0_f64
606  uppert = 1.0_f64 + (rinterval-dx_new*rinterval)
607  call update_v_partial_newton( self, 1.0_f64, dx_new, uppert, lowert, box_new, field_dofs, field )
608  do ind = box_new+1, box_old-1
609  lowert = uppert
610  uppert = lowert + rinterval
611  call update_v_partial_newton( self, 1.0_f64, 0.0_f64, uppert, lowert, ind, field_dofs, field )
612  end do
613  call update_v_partial_newton( self, dx_old, 0.0_f64, 0.0_f64, uppert, box_old, field_dofs, field )
614  field = - field ! Account for the fact that we have integrated the wrong direction
615  end if
616 
617  end subroutine evaluate_int_spline_1d_subnewton
618 
619 
621  subroutine add_current_update_v_spline_pp_1d (self, position_old, position_new, marker_charge, qoverm, bfield_dofs, vi, j_dofs)
622  class(sll_t_particle_mesh_coupling_spline_1d), intent(inout) :: self
623  sll_real64, intent(in) :: position_old(self%dim)
624  sll_real64, intent(in) :: position_new(self%dim)
625  sll_real64, intent(in) :: marker_charge
626  sll_real64, intent(in) :: qoverm
627  sll_real64, intent(in) :: bfield_dofs(:)
628  sll_real64, intent(inout) :: vi(:)
629  sll_real64, intent(inout) :: j_dofs(:)
630  ! local variables
631  sll_int32 :: index_old, index_new, ind
632  sll_real64 :: r_old, r_new
633 
634  call convert_x_to_xbox( self, position_old, r_old, index_old )
635  call convert_x_to_xbox( self, position_new, r_new, index_new )
636 
637  index_old = index_old - 1
638  index_new = index_new - 1
639 
640  if (index_old == index_new) then
641  call update_jv_pp( self, r_old, r_new, index_old, marker_charge, &
642  qoverm, vi(2), j_dofs, bfield_dofs)
643  elseif (index_old < index_new) then
644  call update_jv_pp ( self, r_old, 1.0_f64, index_old, marker_charge, &
645  qoverm, vi(2), j_dofs, bfield_dofs)
646  call update_jv_pp ( self, 0.0_f64, r_new, index_new, marker_charge, &
647  qoverm, vi(2), j_dofs, bfield_dofs)
648  do ind = index_old+1, index_new-1
649  call update_jv_pp ( self, 0.0_f64, 1.0_f64, ind, marker_charge, &
650  qoverm, vi(2), j_dofs, bfield_dofs)
651  end do
652  else
653  call update_jv_pp ( self, 1.0_f64,r_new, index_new, marker_charge, qoverm, &
654  vi(2), j_dofs, bfield_dofs)
655  call update_jv_pp ( self, r_old, 0.0_f64, index_old, marker_charge, qoverm, &
656  vi(2), j_dofs, bfield_dofs)
657  do ind = index_new+1, index_old-1
658  call update_jv_pp ( self, 1.0_f64,0.0_f64, ind, marker_charge, qoverm, &
659  vi(2), j_dofs, bfield_dofs)
660  end do
661  end if
662 
663  end subroutine add_current_update_v_spline_pp_1d
664 
666  subroutine update_jv_pp(self, lower, upper, index, marker_charge, qoverm, vi, j_dofs, bfield_dofs)
667  class(sll_t_particle_mesh_coupling_spline_1d), intent(inout) :: self
668  sll_real64, intent(in) :: lower
669  sll_real64, intent(in) :: upper
670  sll_int32, intent(in) :: index
671  sll_real64, intent(in) :: marker_charge
672  sll_real64, intent(in) :: qoverm
673  sll_real64, intent(inout) :: vi
674  sll_real64, intent(in) :: bfield_dofs(self%n_cells)
675  sll_real64, intent(inout) :: j_dofs(self%n_cells)
676  !Local variables
677  sll_int32 :: ind, i_grid, i_mod, n_cells
678 
679  n_cells = self%n_cells
680  !Evaluation of the primitive integral at the lower and upper bound of the gridcell
681  call sll_s_spline_pp_horner_primitive_1d(self%spline_val, self%spline_degree, self%spline_pp%poly_coeffs_fp, lower)
682  call sll_s_spline_pp_horner_primitive_1d(self%spline_val_more, self%spline_degree, self%spline_pp%poly_coeffs_fp, upper)
683 
684  self%spline_val = (self%spline_val_more - self%spline_val) *(self%delta_x)
685 
686 
687  ind = 1
688  do i_grid = index - self%spline_degree , index
689  i_mod = modulo(i_grid, n_cells ) + 1
690  j_dofs(i_mod) = j_dofs(i_mod) + &
691  (marker_charge*self%spline_val(ind)* self%scaling)
692  vi = vi - qoverm* self%spline_val(ind)*bfield_dofs(i_mod)
693  ind = ind + 1
694  end do
695 
696  end subroutine update_jv_pp
697 
698  !---------------------------------------------------------------------------!
700  subroutine add_current_update_v_spline_1d (self, position_old, position_new, marker_charge, qoverm, bfield_dofs, vi, j_dofs)
701  class(sll_t_particle_mesh_coupling_spline_1d), intent(inout) :: self
702  sll_real64, intent(in) :: position_old(self%dim)
703  sll_real64, intent(in) :: position_new(self%dim)
704  sll_real64, intent(in) :: marker_charge
705  sll_real64, intent(in) :: qoverm
706  sll_real64, intent(in) :: bfield_dofs(:)
707  sll_real64, intent(inout) :: vi(:)
708  sll_real64, intent(inout) :: j_dofs(:)
709  ! local variables
710  sll_int32 :: index_old, index_new, ind
711  sll_real64 :: r_old, r_new
712 
713  call convert_x_to_xbox( self, position_old, r_old, index_old )
714  call convert_x_to_xbox( self, position_new, r_new, index_new )
715 
716  index_old = index_old - 1
717  index_new = index_new - 1
718 
719  if (index_old == index_new) then
720  if (r_old < r_new) then
721  call update_jv( self, r_old, r_new, index_old, marker_charge, &
722  qoverm, 1.0_f64, vi(2), j_dofs, bfield_dofs)
723  else
724  call update_jv( self, r_new, r_old, index_old, marker_charge, qoverm, &
725  -1.0_f64, vi(2), j_dofs, bfield_dofs)
726  end if
727  elseif (index_old < index_new) then
728  call update_jv ( self, r_old, 1.0_f64, index_old, marker_charge, &
729  qoverm, 1.0_f64, vi(2), j_dofs, bfield_dofs)
730  call update_jv ( self, 0.0_f64, r_new, index_new, marker_charge, &
731  qoverm, 1.0_f64, vi(2), j_dofs, bfield_dofs)
732  do ind = index_old+1, index_new-1
733  call update_jv ( self, 0.0_f64, 1.0_f64, ind, marker_charge, &
734  qoverm, 1.0_f64, vi(2), j_dofs, bfield_dofs)
735  end do
736  else
737  call update_jv ( self, r_new, 1.0_f64, index_new, marker_charge, qoverm, &
738  -1.0_f64, vi(2), j_dofs, bfield_dofs)
739  call update_jv ( self, 0.0_f64, r_old, index_old, marker_charge, qoverm, &
740  -1.0_f64, vi(2), j_dofs, bfield_dofs)
741  do ind = index_new+1, index_old-1
742  call update_jv ( self, 0.0_f64, 1.0_f64, ind, marker_charge, qoverm, &
743  -1.0_f64, vi(2), j_dofs, bfield_dofs)
744  end do
745  end if
746 
747 
748  end subroutine add_current_update_v_spline_1d
749 
751  subroutine update_jv(self, lower, upper, index, marker_charge, qoverm, sign, vi, j_dofs, bfield_dofs)
752  class(sll_t_particle_mesh_coupling_spline_1d), intent(inout) :: self
753  sll_real64, intent(in) :: lower
754  sll_real64, intent(in) :: upper
755  sll_int32, intent(in) :: index
756  sll_real64, intent(in) :: marker_charge
757  sll_real64, intent(in) :: qoverm
758  sll_real64, intent(in) :: sign
759  sll_real64, intent(inout) :: vi
760  sll_real64, intent(in) :: bfield_dofs(self%n_cells)
761  sll_real64, intent(inout) :: j_dofs(self%n_cells)
762  !Local variables
763  sll_int32 :: ind, i_grid, i_mod, n_cells, j
764  sll_real64 :: c1, c2
765 
766  n_cells = self%n_cells
767 
768  c1 = 0.5_f64*(upper-lower)
769  c2 = 0.5_f64*(upper+lower)
770 
771 
772  call sll_s_uniform_bsplines_eval_basis(self%spline_degree, c1*self%quad_xw(1,1)+c2, &
773  self%spline_val)
774  self%spline_val = self%spline_val * (self%quad_xw(2,1)*c1)
775  do j=2,self%n_quad_points
776  call sll_s_uniform_bsplines_eval_basis(self%spline_degree, c1*self%quad_xw(1,j)+c2, &
777  self%spline_val_more)
778  self%spline_val = self%spline_val + self%spline_val_more * (self%quad_xw(2,j)*c1)
779  end do
780  self%spline_val = self%spline_val * (sign*self%delta_x)
781 
782  ind = 1
783  do i_grid = index - self%spline_degree , index
784  i_mod = modulo(i_grid, n_cells ) + 1
785  j_dofs(i_mod) = j_dofs(i_mod) + &
786  (marker_charge*self%spline_val(ind)* self%scaling)
787  vi = vi - qoverm* self%spline_val(ind)*bfield_dofs(i_mod)
788  ind = ind + 1
789  end do
790 
791  end subroutine update_jv
792 
794  subroutine update_v_partial(self, upper, lower, uppert, lowert, index, field_dofs, bint )
795  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
796  sll_real64, intent( in ) :: upper
797  sll_real64, intent( in ) :: lower
798  sll_real64, intent( in ) :: uppert
799  sll_real64, intent( in ) :: lowert
800  sll_int32, intent( in ) :: index
801  sll_real64, intent( in ) :: field_dofs(self%n_cells)
802  sll_real64, intent( inout ) :: bint
803  !local variables
804  sll_int32 :: j, ind, i_mod, i_grid, n_cells
805  sll_real64 :: c1, c2, c1t, c2t
806 
807 
808  n_cells = self%n_cells
809 
810  c1 = 0.5_f64*(upper-lower)
811  c2 = 0.5_f64*(upper+lower)
812 
813  c1t = 0.5_f64*(uppert-lowert)
814  c2t = 0.5_f64*(uppert+lowert)
815 
816  call sll_s_uniform_bsplines_eval_basis ( self%spline_degree, c1*self%quadp1_xw(1,1) + c2, &
817  self%spline_val)
818  self%spline_val = self%spline_val * (self%quadp1_xw(2,1)*c1t) * ( c1t* self%quadp1_xw(1,1) + c2t )
819 
820  do j=2, self%n_quadp1_points
821  call sll_s_uniform_bsplines_eval_basis( self%spline_degree, c1*self%quadp1_xw(1,j)+c2, &
822  self%spline_val_more )
823  self%spline_val = self%spline_val + self%spline_val_more * (self%quadp1_xw(2,j)*c1t) * ( c1t* self%quadp1_xw(1,j) + c2t )
824  end do
825  self%spline_val = self%spline_val !* self%delta_x(1)
826 
827  ind = 1
828  do i_grid = index - self%spline_degree, index
829  i_mod = modulo(i_grid-1, n_cells ) + 1
830  bint = bint + self%spline_val(ind) * field_dofs( i_mod )
831  ind = ind+1
832  end do
833 
834  end subroutine update_v_partial
835 
836 
837 
839  subroutine update_v_partial_newton(self, upper, lower, uppert, lowert, index, field_dofs, bint )
840  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
841  sll_real64, intent( in ) :: upper
842  sll_real64, intent( in ) :: lower
843  sll_real64, intent( in ) :: uppert
844  sll_real64, intent( in ) :: lowert
845  sll_int32, intent( in ) :: index
846  sll_real64, intent( in ) :: field_dofs(self%n_cells)
847  sll_real64, intent( inout ) :: bint
848 
849  sll_int32 :: j, ind, i_mod, i_grid, n_cells
850  sll_real64 :: c1, c2, c1t, c2t, tau
851 
852  n_cells = self%n_cells
853 
854  c1 = 0.5_f64*(upper-lower)
855  c2 = 0.5_f64*(upper+lower)
856 
857  c1t = 0.5_f64*(uppert-lowert)
858  c2t = 0.5_f64*(uppert+lowert)
859 
860 
861  call sll_s_uniform_bsplines_eval_deriv ( self%spline_degree, c1*self%quadp1_xw(1,1) + c2, &
862  self%spline_val)
863  tau = ( c1t* self%quadp1_xw(1,1) + c2t )
864  self%spline_val = self%spline_val * (self%quadp1_xw(2,1)*c1t) * ( 1.0_f64 - tau ) * tau
865 
866  do j=2, self%n_quadp1_points
867  call sll_s_uniform_bsplines_eval_deriv( self%spline_degree, c1*self%quadp1_xw(1,j)+c2, &
868  self%spline_val_more )
869  tau = ( c1t* self%quadp1_xw(1,j) + c2t )
870  self%spline_val = self%spline_val + self%spline_val_more * (self%quadp1_xw(2,j)*c1t) * ( 1.0_f64 - tau ) * tau
871  end do
872  self%spline_val = self%spline_val !* self%delta_x(1)
873 
874  ind = 1
875  do i_grid = index - self%spline_degree, index
876  i_mod = modulo(i_grid-1, n_cells ) + 1
877  bint = bint + self%spline_val(ind) * field_dofs( i_mod )
878  ind = ind+1
879  end do
880 
881  end subroutine update_v_partial_newton
882 
883 
885  subroutine convert_x_to_xbox( self, position, xi, box )
886  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
887  sll_real64, intent( in ) :: position(self%dim)
888  sll_real64, intent( out ) :: xi
889  sll_int32, intent( out ) :: box
890 
891  xi = (position(1) - self%domain(1)) /self%delta_x
892  box = floor( xi ) + 1
893  xi = xi - real(box-1, f64)
894 
895  end subroutine convert_x_to_xbox
896 
897 
899  subroutine init_spline_1d( self, domain, n_cells, no_particles, spline_degree, smoothing_type )
900  class( sll_t_particle_mesh_coupling_spline_1d), intent(out) :: self
901  sll_int32, intent(in) :: n_cells
902  sll_real64, intent(in) :: domain(2)
903  sll_int32, intent(in) :: no_particles
904  sll_int32, intent(in) :: spline_degree
905  sll_int32, intent(in) :: smoothing_type
906  !local variables
907  sll_int32 :: ierr
908 
909 
910  self%dim = 1
911 
912  ! Store grid information
913  self%domain = domain
914  self%n_cells = n_cells
915  self%n_dofs = n_cells
916  self%delta_x = (self%domain(2)-self%domain(1))/real(n_cells, f64)
917 
918  ! Store basis function information
919  self%no_particles = no_particles
920 
921  ! Initialize information on the spline
922  self%spline_degree = spline_degree
923  self%n_span = spline_degree + 1
924 
925  ! Initialize information on smoothing type
926  if (smoothing_type == sll_p_collocation) then
927  self%scaling = 1.0_f64/self%delta_x
928  elseif (smoothing_type == sll_p_galerkin) then
929  self%scaling = 1.0_f64
930  else
931  print*, 'Smoothing Type ', smoothing_type, ' not implemented for kernel_smoother_spline_1d.'
932  end if
933 
934  self%n_quad_points = (self%spline_degree+2)/2
935 
936  ALLOCATE( self%spline_val(self%n_span), stat = ierr)
937  sll_assert( ierr == 0 )
938  ALLOCATE( self%spline_val_more(self%n_span), stat = ierr )
939  sll_assert( ierr == 0 )
940  ALLOCATE( self%quad_xw(2,self%n_quad_points), stat = ierr )
941  sll_assert( ierr == 0 )
942 
943  ! normalized Gauss Legendre points and weights
944  self%quad_xw = sll_f_gauss_legendre_points_and_weights(self%n_quad_points)
945  call sll_s_spline_pp_init_1d( self%spline_pp, spline_degree, n_cells)
946 
947  ! For nonlinear disgradE
948  self%n_quadp1_points = (self%spline_degree+3)/2
949  allocate( self%quadp1_xw(2, self%n_quadp1_points) )
950  allocate( self%spline_pol_val(self%n_span) )
951  ! normalized Gauss Legendre points and weights
952  self%quadp1_xw = sll_f_gauss_legendre_points_and_weights(self%n_quadp1_points)
953 
954  ! For implicit subcycling
955  self%n_quadp2_points = (self%spline_degree+4)/2
956  allocate( self%quadp2_xw(2, self%n_quadp2_points) )
957  ! normalized Gauss Legendre points and weights
958  self%quadp2_xw = sll_f_gauss_legendre_points_and_weights(self%n_quadp2_points)
959 
960 
961  ! For smoothed evaluate and add_charge
962  allocate( self%quads1_xw(2, self%n_quadp1_points) )
963  ! normalized Gauss Legendre points and weights to [0,1]
964  self%quads1_xw = self%quadp1_xw
965  self%quads1_xw(1,:) = (self%quads1_xw(1,:) + 1.0_f64)*0.5_f64
966  self%quads1_xw(2,:) = self%quads1_xw(2,:)*0.5_f64
967 
968 
969 
970  ! For smoothed add_current
971  self%n_quads2_points = (self%spline_degree+4)/2
972  allocate( self%quads2_xw(2, self%n_quads2_points) )
973  ! normalized Gauss Legendre points and weights to [0,1]
974  self%quads2_xw = sll_f_gauss_legendre_points_and_weights(self%n_quads2_points)
975  self%quads2_xw(1,:) = (self%quads2_xw(1,:) + 1.0_f64)*0.5_f64
976  self%quads2_xw(2,:) = self%quads2_xw(2,:)*0.5_f64
977 
978  self%n_quad_points_line = (self%spline_degree+2)/2
979 
980  allocate( self%quad_xw_line(2,self%n_quad_points_line) )
981  ! normalized Gauss Legendre points and weights
982  self%quad_xw_line = sll_f_gauss_legendre_points_and_weights(self%n_quad_points_line)
983 
984  self%quad_xw_line(1,:) = 0.5_f64*(self%quad_xw_line(1,:)+1.0_f64)
985  self%quad_xw_line(2,:) = 0.5_f64*self%quad_xw_line(2,:)
986 
987  end subroutine init_spline_1d
988 
990  subroutine free_spline_1d(self)
991  class(sll_t_particle_mesh_coupling_spline_1d), intent( inout ) :: self
992 
993  deallocate(self%spline_val)
994  deallocate(self%spline_val_more)
995  deallocate(self%quad_xw)
996  call sll_s_spline_pp_free_1d(self%spline_pp)
997 
998  end subroutine free_spline_1d
999 
1000 
1001  !-------------------------------------------------------------------------------------------
1002 
1003  subroutine sll_s_new_particle_mesh_coupling_spline_1d_ptr(smoother, domain, n_cells, no_particles, spline_degree, smoothing_type)
1004  class( sll_c_particle_mesh_coupling_1d), pointer, intent(out) :: smoother
1005  sll_int32, intent(in) :: n_cells
1006  sll_real64, intent(in) :: domain(2)
1007  sll_int32, intent(in) :: no_particles
1008  sll_int32, intent(in) :: spline_degree
1009  sll_int32, intent(in) :: smoothing_type
1010 
1011  !local variables
1012  sll_int32 :: ierr
1013 
1014 
1015  sll_allocate( sll_t_particle_mesh_coupling_spline_1d :: smoother , ierr)
1016  sll_assert( ierr == 0)
1017 
1018  select type( smoother )
1020  call smoother%init( domain, n_cells, no_particles, spline_degree, smoothing_type )
1021  end select
1022 
1024 
1025  !----------------------------------------------------------------------------------
1026 
1027  subroutine sll_s_new_particle_mesh_coupling_spline_1d(smoother, domain, n_cells, no_particles, spline_degree, smoothing_type)
1028  class( sll_c_particle_mesh_coupling_1d), allocatable, intent(out):: smoother
1029  sll_int32, intent(in) :: n_cells
1030  sll_real64, intent(in) :: domain(2)
1031  sll_int32, intent(in) :: no_particles
1032  sll_int32, intent(in) :: spline_degree
1033  sll_int32, intent(in) :: smoothing_type
1034 
1035  !local variables
1036  sll_int32 :: ierr
1037 
1038 
1039  allocate( sll_t_particle_mesh_coupling_spline_1d :: smoother , stat=ierr)
1040  sll_assert( ierr == 0)
1041 
1042  select type( smoother )
1044  call smoother%init( domain, n_cells, no_particles, spline_degree, smoothing_type )
1045  end select
1046 
1048 
1049 
1050 
1051 
1053  subroutine add_current_split_spline_1d(self, position_old, position_new, iter, total_iter, marker_charge, j_dofs )
1054  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
1055  sll_real64, intent( in ) :: position_old(self%dim)
1056  sll_real64, intent( in ) :: position_new(self%dim)
1057  sll_int32, intent( in ) :: iter
1058  sll_int32, intent( in ) :: total_iter
1059  sll_real64, intent( in ) :: marker_charge
1060  sll_real64, intent( inout ) :: j_dofs(self%n_cells,2)
1061  !local variables
1062  sll_real64 :: dx_new, dx_old
1063  sll_int32 :: box_new, box_old, ind
1064  sll_real64 :: rinterval, lowert, uppert
1065 
1066  dx_new = (position_new(1)-self%domain(1))/self%delta_x
1067  box_new = floor(dx_new)+1
1068  dx_new = dx_new - real(box_new-1,f64)
1069  dx_old = (position_old(1)-self%domain(1))/self%delta_x
1070  box_old = floor(dx_old)+1
1071  dx_old = dx_old - real(box_old-1,f64)
1072 
1073  if ( box_old == box_new ) then
1074  call update_current_partial(self, dx_new, dx_old, 1.0_f64, 0.0_f64, box_old, &
1075  iter, total_iter, marker_charge, j_dofs )
1076  elseif ( box_old < box_new ) then
1077  rinterval = self%delta_x/(position_new(1)-position_old(1))
1078  lowert = 0.0_f64
1079  uppert = (1.0_f64-dx_old)*rinterval
1080  call update_current_partial( self, 1.0_f64, dx_old, uppert, lowert, box_old, &
1081  iter, total_iter, marker_charge, j_dofs )
1082  do ind = box_old+1, box_new-1
1083  lowert = uppert
1084  uppert = uppert + rinterval
1085  call update_current_partial( self, 1.0_f64, 0.0_f64, uppert, lowert, ind, &
1086  iter, total_iter, marker_charge, j_dofs )
1087  end do
1088  call update_current_partial( self, dx_new, 0.0_f64, 1.0_f64, uppert, box_new, &
1089  iter, total_iter, marker_charge, j_dofs )
1090  else
1091  rinterval = self%delta_x/(position_new(1)-position_old(1))
1092  lowert = 1.0_f64
1093  uppert = 1.0_f64 + (rinterval - dx_new*rinterval)
1094  ! Note the negative sign in the marker_charge accounts for the fact that we are integrating in the wrong direction
1095  call update_current_partial( self, 1.0_f64, dx_new, uppert, lowert, box_new, &
1096  iter, total_iter, -marker_charge, j_dofs )
1097  do ind = box_new+1, box_old-1
1098  lowert = uppert
1099  uppert = lowert + rinterval
1100  call update_current_partial( self, 1.0_f64, 0.0_f64, uppert, lowert, ind, &
1101  iter, total_iter, -marker_charge, j_dofs )
1102  end do
1103  call update_current_partial( self, dx_old, 0.0_f64, 0.0_f64, uppert, box_old, &
1104  iter, total_iter, -marker_charge, j_dofs )
1105  end if
1106 
1107  end subroutine add_current_split_spline_1d
1108 
1109 
1111  subroutine update_current_partial(self, upper, lower, uppert, lowert, index, iter, total_iter, marker_charge, j_dofs )
1112  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
1113  sll_real64, intent( in ) :: upper
1114  sll_real64, intent( in ) :: lower
1115  sll_real64, intent( in ) :: uppert
1116  sll_real64, intent( in ) :: lowert
1117  sll_int32, intent( in ) :: index
1118  sll_int32, intent( in ) :: iter
1119  sll_int32, intent( in ) :: total_iter
1120  sll_real64, intent( in ) :: marker_charge
1121  sll_real64, intent( inout ) :: j_dofs(self%n_cells,2)
1122 
1123 
1124  sll_int32 :: j, ind, i_mod, i_grid, n_cells
1125  sll_real64 :: c1, c2, c1t, c2t
1126  sll_real64 :: time
1127 
1128  n_cells = self%n_cells
1129 
1130  c1 = 0.5_f64*(upper-lower)
1131  c2 = 0.5_f64*(upper+lower)
1132 
1133  c1t = 0.5_f64*(uppert-lowert)
1134  c2t = 0.5_f64*(uppert+lowert)
1135 
1136  call sll_s_uniform_bsplines_eval_basis ( self%spline_degree, &
1137  c1*self%quadp1_xw(1,1) + c2, &
1138  self%spline_val)
1139  time = (c1t* self%quadp1_xw(1,1) + c2t + real(iter,f64))/real(total_iter, f64)
1140  self%spline_pol_val = self%spline_val * (self%quadp1_xw(2,1)*c1t) * &
1141  ( 1.0_f64 - time )
1142  self%spline_val = self%spline_val * (self%quadp1_xw(2,1)*c1t) * time
1143 
1144  do j=2, self%n_quadp1_points
1145  call sll_s_uniform_bsplines_eval_basis( self%spline_degree, &
1146  c1*self%quadp1_xw(1,j)+c2, &
1147  self%spline_val_more )
1148  time = (c1t* self%quadp1_xw(1,j) + c2t + real(iter,f64))/real(total_iter, f64)
1149  self%spline_val = self%spline_val + &
1150  self%spline_val_more * (self%quadp1_xw(2,j)*c1t) * time
1151  self%spline_pol_val = self%spline_pol_val + &
1152  self%spline_val_more * (self%quadp1_xw(2,j)*c1t) * ( 1.0_f64 - time )
1153  end do
1154 
1155  ind = 1
1156  do i_grid = index - self%spline_degree, index
1157  i_mod = modulo(i_grid-1, n_cells ) + 1
1158  j_dofs(i_mod,1) = j_dofs(i_mod,1) + marker_charge * self%spline_val(ind)
1159  j_dofs(i_mod,2) = j_dofs(i_mod,2) + marker_charge * self%spline_pol_val(ind)
1160  ind = ind+1
1161  end do
1162 
1163  end subroutine update_current_partial
1164 
1166  subroutine evaluate_int_linear_quad_spline_1d(self, position_old, position_new, iter, total_iter, field_dofs_1, field_dofs_2, field)
1167  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
1168  sll_real64, intent( in ) :: position_old(self%dim)
1169  sll_real64, intent( in ) :: position_new(self%dim)
1170  sll_real64, intent( in ) :: field_dofs_1(self%n_cells)
1171  sll_real64, intent( in ) :: field_dofs_2(self%n_cells)
1172  sll_int32, intent( in ) :: iter
1173  sll_int32, intent( in ) :: total_iter
1174  sll_real64, intent( out ) :: field(2)
1175  !local variables
1176  sll_real64 :: dx_new, dx_old
1177  sll_int32 :: box_new, box_old, ind
1178  sll_real64 :: rinterval, lowert, uppert
1179 
1180  dx_new = (position_new(1)-self%domain(1))/self%delta_x
1181  box_new = floor(dx_new)+1
1182  dx_new = dx_new - real(box_new-1,f64)
1183  dx_old = (position_old(1)-self%domain(1))/self%delta_x
1184  box_old = floor(dx_old)+1
1185  dx_old = dx_old - real(box_old-1,f64)
1186 
1187  if ( box_old == box_new ) then
1188  rinterval = 0.0_f64
1189  else
1190  rinterval = self%delta_x/(position_new(1)-position_old(1))
1191  end if
1192  field = 0.0_f64
1193 
1194  if ( box_old == box_new ) then
1195  call update_field_linear_partial(self, dx_new, dx_old, 1.0_f64, 0.0_f64, box_old, &
1196  field_dofs_1, field_dofs_2, iter, total_iter, field )
1197  elseif ( box_old < box_new ) then
1198  lowert = 0.0_f64
1199  uppert = (1.0_f64-dx_old)*rinterval
1200  call update_field_linear_partial( self, 1.0_f64, dx_old, uppert, lowert, box_old, &
1201  field_dofs_1, field_dofs_2, iter, total_iter, field )
1202  do ind = box_old+1, box_new-1
1203  lowert = uppert
1204  uppert = uppert + rinterval
1205  call update_field_linear_partial( self, 1.0_f64, 0.0_f64, uppert, lowert, ind, &
1206  field_dofs_1, field_dofs_2, iter, total_iter, field )
1207  end do
1208  call update_field_linear_partial( self, dx_new, 0.0_f64, 1.0_f64, uppert, box_new, &
1209  field_dofs_1, field_dofs_2, iter, total_iter, field )
1210  else
1211  lowert = 1.0_f64
1212  uppert = 1.0_f64 + (rinterval-dx_new*rinterval)
1213  call update_field_linear_partial( self, 1.0_f64, dx_new, uppert, lowert, box_new, &
1214  field_dofs_1, field_dofs_2, iter, total_iter, field )
1215  do ind = box_new+1, box_old-1
1216  lowert = uppert
1217  uppert = lowert + rinterval
1218  call update_field_linear_partial( self, 1.0_f64, 0.0_f64, uppert, lowert, ind, &
1219  field_dofs_1, field_dofs_2, iter, total_iter, field )
1220  end do
1221  call update_field_linear_partial( self, dx_old, 0.0_f64, 0.0_f64, uppert, box_old, &
1222  field_dofs_1, field_dofs_2, iter, total_iter, field )
1223  field = - field ! Account for the fact that we have integrated the wrong direction
1224  end if
1225 
1226  end subroutine evaluate_int_linear_quad_spline_1d
1227 
1228 
1230  subroutine update_field_linear_partial(self, upper, lower, uppert, lowert, index, field_dofs_1, field_dofs_2, iter, total_iter, bint )
1231  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
1232  sll_real64, intent( in ) :: upper
1233  sll_real64, intent( in ) :: lower
1234  sll_real64, intent( in ) :: uppert
1235  sll_real64, intent( in ) :: lowert
1236  sll_int32, intent( in ) :: index
1237  sll_real64, intent( in ) :: field_dofs_1(self%n_cells)
1238  sll_real64, intent( in ) :: field_dofs_2(self%n_cells)
1239  sll_int32, intent( in ) :: iter
1240  sll_int32, intent( in ) :: total_iter
1241  sll_real64, intent( inout ) :: bint(2)
1242 
1243  sll_int32 :: j, ind, i_mod, i_grid, n_cells
1244  sll_real64 :: c1, c2, c1t, c2t
1245  sll_real64 :: tau, time, weight, field_dof_t
1246 
1247  n_cells = self%n_cells
1248 
1249  c1 = 0.5_f64*(upper-lower)
1250  c2 = 0.5_f64*(upper+lower)
1251 
1252  c1t = 0.5_f64*(uppert-lowert)
1253  c2t = 0.5_f64*(uppert+lowert)
1254 
1255  do j=1, self%n_quadp2_points
1256  tau = c1t* self%quadp2_xw(1,j) + c2t
1257  time = ( tau +real(iter,f64))/real(total_iter, f64)
1258  weight = (self%quadp2_xw(2,j)*c1t)
1259  call sll_s_uniform_bsplines_eval_basis( self%spline_degree, c1*self%quadp2_xw(1,j)+c2, &
1260  self%spline_val )
1261  ind = 1
1262  do i_grid = index - self%spline_degree, index
1263  i_mod = modulo(i_grid-1, n_cells ) + 1
1264  field_dof_t = field_dofs_1( i_mod ) + time * ( field_dofs_2( i_mod ) - field_dofs_1( i_mod ) )
1265  bint(1) = bint(1) + self%spline_val(ind) * field_dof_t * weight * tau
1266  bint(2) = bint(2) + self%spline_val(ind) * field_dof_t * weight * (1.0_f64 - tau )
1267  ind = ind+1
1268  end do
1269  end do
1270 
1271  end subroutine update_field_linear_partial
1272 
1274  subroutine evaluate_int_linear_quad_subnewton_spline_1d(self, position_old, position_new, iter, total_iter, field_dofs_1, field_dofs_2, field)
1275  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
1276  sll_real64, intent( in ) :: position_old(self%dim)
1277  sll_real64, intent( in ) :: position_new(self%dim)
1278  sll_real64, intent( in ) :: field_dofs_1(self%n_cells)
1279  sll_real64, intent( in ) :: field_dofs_2(self%n_cells)
1280  sll_int32, intent( in ) :: iter
1281  sll_int32, intent( in ) :: total_iter
1282  sll_real64, intent( out ) :: field
1283  !local variables
1284  sll_real64 :: dx_new, dx_old
1285  sll_int32 :: box_new, box_old, ind
1286  sll_real64 :: rinterval, lowert, uppert
1287 
1288  dx_new = (position_new(1)-self%domain(1))/self%delta_x
1289  box_new = floor(dx_new)+1
1290  dx_new = dx_new - real(box_new-1,f64)
1291  dx_old = (position_old(1)-self%domain(1))/self%delta_x
1292  box_old = floor(dx_old)+1
1293  dx_old = dx_old - real(box_old-1,f64)
1294 
1295  rinterval = self%delta_x/(position_new(1)-position_old(1))
1296 
1297  field = 0.0_f64
1298 
1299  if ( box_old == box_new ) then
1300  call update_field_linear_partial_newton(self, dx_new, dx_old, 1.0_f64, 0.0_f64, box_old, &
1301  field_dofs_1, field_dofs_2, iter, total_iter, field )
1302  elseif ( box_old < box_new ) then
1303  lowert = 0.0_f64
1304  uppert = (1.0_f64-dx_old)*rinterval
1305  call update_field_linear_partial_newton( self, 1.0_f64, dx_old, uppert, lowert, box_old, &
1306  field_dofs_1, field_dofs_2, iter, total_iter, field )
1307  do ind = box_old+1, box_new-1
1308  lowert = uppert
1309  uppert = uppert + rinterval
1310  call update_field_linear_partial_newton( self, 1.0_f64, 0.0_f64, uppert, lowert, ind, &
1311  field_dofs_1, field_dofs_2, iter, total_iter, field )
1312  end do
1313  call update_field_linear_partial_newton( self, dx_new, 0.0_f64, 1.0_f64, uppert, box_new, &
1314  field_dofs_1, field_dofs_2, iter, total_iter, field )
1315  else
1316  lowert = 1.0_f64
1317  uppert = 1.0_f64 + (rinterval-dx_new*rinterval)
1318  call update_field_linear_partial_newton( self, 1.0_f64, dx_new, uppert, lowert, box_new, &
1319  field_dofs_1, field_dofs_2, iter, total_iter, field )
1320  do ind = box_new+1, box_old-1
1321  lowert = uppert
1322  uppert = lowert + rinterval
1323  call update_field_linear_partial_newton( self, 1.0_f64, 0.0_f64, uppert, lowert, ind, &
1324  field_dofs_1, field_dofs_2, iter, total_iter, field )
1325  end do
1326  call update_field_linear_partial_newton( self, dx_old, 0.0_f64, 0.0_f64, uppert, box_old, &
1327  field_dofs_1, field_dofs_2, iter, total_iter, field )
1328  field = - field ! Account for the fact that we have integrated the wrong direction
1329  end if
1330 
1332 
1333 
1335  subroutine update_field_linear_partial_newton(self, upper, lower, uppert, lowert, index, field_dofs_1, field_dofs_2, iter, total_iter, bint )
1336  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
1337  sll_real64, intent( in ) :: upper
1338  sll_real64, intent( in ) :: lower
1339  sll_real64, intent( in ) :: uppert
1340  sll_real64, intent( in ) :: lowert
1341  sll_int32, intent( in ) :: index
1342  sll_real64, intent( in ) :: field_dofs_1(self%n_cells)
1343  sll_real64, intent( in ) :: field_dofs_2(self%n_cells)
1344  sll_int32, intent( in ) :: iter
1345  sll_int32, intent( in ) :: total_iter
1346  sll_real64, intent( inout ) :: bint
1347 
1348  sll_int32 :: j, ind, i_mod, i_grid, n_cells
1349  sll_real64 :: c1, c2, c1t, c2t
1350  sll_real64 :: tau, time, weight, field_dof_t
1351 
1352  n_cells = self%n_cells
1353 
1354  c1 = 0.5_f64*(upper-lower)
1355  c2 = 0.5_f64*(upper+lower)
1356 
1357  c1t = 0.5_f64*(uppert-lowert)
1358  c2t = 0.5_f64*(uppert+lowert)
1359 
1360  do j=1, self%n_quadp2_points
1361  tau = c1t* self%quadp2_xw(1,j) + c2t
1362  time = ( tau +real(iter,f64))/real(total_iter, f64)
1363  weight = (self%quadp2_xw(2,j)*c1t)
1364  call sll_s_uniform_bsplines_eval_deriv ( self%spline_degree, c1*self%quadp2_xw(1,1) + c2, &
1365  self%spline_val)
1366  ind = 1
1367  do i_grid = index - self%spline_degree, index
1368  i_mod = modulo(i_grid-1, n_cells ) + 1
1369  field_dof_t = field_dofs_1( i_mod ) + time * ( field_dofs_2( i_mod ) - field_dofs_1( i_mod ) )
1370  bint = bint + self%spline_val(ind) * field_dof_t * weight * (1.0_f64 - tau ) * tau
1371  ind = ind+1
1372  end do
1373  end do
1374 
1375  end subroutine update_field_linear_partial_newton
1376 
1378  subroutine add_charge_int_spline_1d(self, position_old, position_new, marker_charge, j_dofs )
1379  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
1380  sll_real64, intent( in ) :: position_old(self%dim)
1381  sll_real64, intent( in ) :: position_new(self%dim)
1382  sll_real64, intent( in ) :: marker_charge
1383  sll_real64, intent( inout ) :: j_dofs(self%n_cells)
1384  !local variables
1385  sll_real64 :: dx_new, dx_old
1386  sll_int32 :: box_new, box_old, ind
1387  sll_real64 :: rinterval, lowert, uppert
1388 
1389  dx_new = (position_new(1)-self%domain(1))/self%delta_x
1390  box_new = floor(dx_new)+1
1391  dx_new = dx_new - real(box_new-1,f64)
1392  dx_old = (position_old(1)-self%domain(1))/self%delta_x
1393  box_old = floor(dx_old)+1
1394  dx_old = dx_old - real(box_old-1,f64)
1395 
1396  if ( box_old == box_new ) then
1397  call update_charge_int_partial(self, dx_new, dx_old, 1.0_f64, 0.0_f64, box_old, &
1398  marker_charge, j_dofs )
1399  elseif ( box_old < box_new ) then
1400  rinterval = self%delta_x/(position_new(1)-position_old(1))
1401  lowert = 0.0_f64
1402  uppert = (1.0_f64-dx_old)*rinterval
1403  call update_charge_int_partial( self, 1.0_f64, dx_old, uppert, lowert, box_old, &
1404  marker_charge, j_dofs )
1405  do ind = box_old+1, box_new-1
1406  lowert = uppert
1407  uppert = uppert + rinterval
1408  call update_charge_int_partial( self, 1.0_f64, 0.0_f64, uppert, lowert, ind, &
1409  marker_charge, j_dofs )
1410  end do
1411  call update_charge_int_partial( self, dx_new, 0.0_f64, 1.0_f64, uppert, box_new, &
1412  marker_charge, j_dofs )
1413  else
1414  rinterval = self%delta_x/(position_new(1)-position_old(1))
1415  lowert = 1.0_f64
1416  uppert = 1.0_f64 + (rinterval - dx_new*rinterval)
1417  ! Note the negative sign in the marker_charge accounts for the fact that we are integrating in the wrong direction
1418  call update_charge_int_partial( self, 1.0_f64, dx_new, uppert, lowert, box_new, &
1419  -marker_charge, j_dofs )
1420  do ind = box_new+1, box_old-1
1421  lowert = uppert
1422  uppert = lowert + rinterval
1423  call update_charge_int_partial( self, 1.0_f64, 0.0_f64, uppert, lowert, ind, &
1424  -marker_charge, j_dofs )
1425  end do
1426  call update_charge_int_partial( self, dx_old, 0.0_f64, 0.0_f64, uppert, box_old, &
1427  -marker_charge, j_dofs )
1428  end if
1429 
1430  end subroutine add_charge_int_spline_1d
1431 
1432 
1434  subroutine update_charge_int_partial(self, upper, lower, uppert, lowert, index, marker_charge, j_dofs )
1435  class( sll_t_particle_mesh_coupling_spline_1d ), intent(inout) :: self
1436  sll_real64, intent( in ) :: upper
1437  sll_real64, intent( in ) :: lower
1438  sll_real64, intent( in ) :: uppert
1439  sll_real64, intent( in ) :: lowert
1440  sll_int32, intent( in ) :: index
1441  sll_real64, intent( in ) :: marker_charge
1442  sll_real64, intent( inout ) :: j_dofs(self%n_cells)
1443  !local variables
1444  sll_int32 :: j, ind, i_mod, i_grid, n_cells
1445  sll_real64 :: c1, c2, c1t, c2t
1446 
1447  n_cells = self%n_cells
1448 
1449  c1 = 0.5_f64*(upper-lower)
1450  c2 = 0.5_f64*(upper+lower)
1451 
1452  c1t = 0.5_f64*(uppert-lowert)
1453  c2t = 0.5_f64*(uppert+lowert)
1454 
1455  call sll_s_uniform_bsplines_eval_basis ( self%spline_degree, &
1456  c1*self%quadp1_xw(1,1) + c2, &
1457  self%spline_val)
1458  self%spline_val = self%spline_val * (self%quadp1_xw(2,1)*c1t)
1459 
1460  do j=2, self%n_quadp1_points
1461  call sll_s_uniform_bsplines_eval_basis( self%spline_degree, &
1462  c1*self%quadp1_xw(1,j)+c2, &
1463  self%spline_val_more )
1464  self%spline_val = self%spline_val + &
1465  self%spline_val_more * (self%quadp1_xw(2,j)*c1t)
1466  end do
1467 
1468  ind = 1
1469  do i_grid = index - self%spline_degree, index
1470  i_mod = modulo(i_grid-1, n_cells ) + 1
1471  j_dofs(i_mod) = j_dofs(i_mod) + marker_charge * self%spline_val(ind)
1472  ind = ind+1
1473  end do
1474 
1475  end subroutine update_charge_int_partial
1476 
1477 
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 [...
Low level arbitrary degree splines.
Base class for kernel smoothers for accumulation and field evaluation in PIC.
integer(kind=i32), parameter, public sll_p_collocation
Kernel smoother for 1d with splines of arbitrary degree placed on a uniform mesh.
subroutine add_particle_mass_spline_1d(self, position, marker_charge, particle_mass)
Add the contribution of one particle to the approximate mass matrix.
subroutine update_charge_int_partial(self, upper, lower, uppert, lowert, index, marker_charge, j_dofs)
Helper function for add_charge_int (takes care of the per cell computations)
subroutine update_field_linear_partial_newton(self, upper, lower, uppert, lowert, index, field_dofs_1, field_dofs_2, iter, total_iter, bint)
Helper function for evaluate_int_linear_quad_subnewton (takes care of the per cell computations)
subroutine add_particle_mass_spline_1d_full(self, position, marker_charge, particle_mass)
Add the contribution of one particle to the approximate mass matrix.
subroutine add_charge_single_spline_1d(self, position, marker_charge, rho_dofs)
helper function to compute the integral of j using Gauss quadrature (needs to be provided for the der...
subroutine add_charge_int_spline_1d(self, position_old, position_new, marker_charge, j_dofs)
For explicit version of the subcycling propagator (to evaluate the charge as it comes out integrated ...
subroutine evaluate_int_quad_spline_1d(self, position_old, position_new, field_dofs, field)
Evaluates an integrated field (between position_old and position_new) but with a quadrature formula f...
subroutine update_v_partial_newton(self, upper, lower, uppert, lowert, index, field_dofs, bint)
Helper function for evaluate_int_quad_subnewton (despite the name it evaluates a field and does not u...
subroutine add_current_spline_1d(self, position_old, position_new, marker_charge, j_dofs)
Add current with integration over x.
subroutine evaluate_field_single_spline_1d(self, position, field_dofs, field_value)
Evaluate field at at position position.
subroutine update_jv(self, lower, upper, index, marker_charge, qoverm, sign, vi, j_dofs, bfield_dofs)
Helper function for add_current_update_v.
subroutine add_current_update_v_spline_pp_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)
subroutine, public sll_s_new_particle_mesh_coupling_spline_1d_ptr(smoother, domain, n_cells, no_particles, spline_degree, smoothing_type)
subroutine update_jv_pp(self, lower, upper, index, marker_charge, qoverm, vi, j_dofs, bfield_dofs)
Helper function for add_current_update_v.
subroutine add_current_update_v_spline_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)
subroutine init_spline_1d(self, domain, n_cells, no_particles, spline_degree, smoothing_type)
Initializer.
subroutine evaluate_multiple_spline_1d(self, position, components, field_dofs, field_value)
Evaluate several fields at position position.
subroutine evaluate_int_spline_1d_subnewton(self, position_old, position_new, field_dofs, field)
subroutine evaluate_int_spline_1d(self, position_old, position_new, vbar, field_dofs, field)
Evaluates the integral int_{poisition_old}^{position_new} field(x) d x.
subroutine add_current_evaluate_spline_1d(self, position_old, position_new, marker_charge, vbar, field_dofs, j_dofs, field)
Combines add_current and evaluate_int.
subroutine, public sll_s_new_particle_mesh_coupling_spline_1d(smoother, domain, n_cells, no_particles, spline_degree, smoothing_type)
subroutine convert_x_to_xbox(self, position, xi, box)
Helper function that identifies in which box the particle is found and its normalized poistion in the...
subroutine update_current_partial(self, upper, lower, uppert, lowert, index, iter, total_iter, marker_charge, j_dofs)
Helper function for add_current_split.
subroutine evaluate_int_linear_quad_spline_1d(self, position_old, position_new, iter, total_iter, field_dofs_1, field_dofs_2, field)
Evaluates the integrals with tau function for implicit subcycling.
subroutine add_current_split_spline_1d(self, position_old, position_new, iter, total_iter, marker_charge, j_dofs)
Add current version with tau and (1-tau) which is needed for the implicit subcycling scheme.
subroutine update_v_partial(self, upper, lower, uppert, lowert, index, field_dofs, bint)
Helper function for evaluate_int_quad (despite the name it evaluates a field and does not update v)
subroutine evaluate_int_linear_quad_subnewton_spline_1d(self, position_old, position_new, iter, total_iter, field_dofs_1, field_dofs_2, field)
For implicit subcycling propagator, evaluates the integral needed for Newtons's method in the evaluat...
subroutine update_field_linear_partial(self, upper, lower, uppert, lowert, index, field_dofs_1, field_dofs_2, iter, total_iter, bint)
Helper function for evaluate_int_linear_quad (implements part within one cell)
Splines in pp form.
subroutine, public sll_s_spline_pp_horner_primitive_1d(val, degree, pp_coeffs, x)
Perform a 1d hornerschema on the pp_coeffs evaluate at x.
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)
subroutine, public sll_s_spline_pp_free_1d(self)
Destructor 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.
    Report Typos and Errors