sll_m_particle_mesh_coupling_spline_cl_1d.F90

Go to the documentation of this file.

 
module sll_m_particle_mesh_coupling_spline_cl_1d
 
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_assert.h"
#include "sll_memory.h"
#include "sll_working_precision.h"
#include "sll_errors.h"
 
  use sll_m_low_level_bsplines, only: &
       sll_s_uniform_bsplines_eval_basis
 
  use sll_m_gauss_legendre_integration, only : &
       sll_f_gauss_legendre_points_and_weights
 
  use sll_m_particle_mesh_coupling_base_1d, only: &
       sll_p_collocation, &
       sll_c_particle_mesh_coupling_1d, &
       sll_p_galerkin
 
  use sll_m_particle_group_base, only: &
       sll_c_particle_group_base
 
  use sll_m_splines_pp, only: &
       sll_t_spline_pp_1d, &
       sll_s_spline_pp_horner_m_1d, &
       sll_s_spline_pp_horner_primitive_1d, &
       sll_f_spline_pp_horner_1d, &
       sll_s_spline_pp_free_1d, &
       sll_s_spline_pp_init_1d
  implicit none
 
  public :: &
       sll_t_particle_mesh_coupling_spline_cl_1d, &
       sll_s_new_particle_mesh_coupling_spline_cl_1d_ptr, &
       sll_s_new_particle_mesh_coupling_spline_cl_1d
 
  private
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  type, extends(sll_c_particle_mesh_coupling_1d) :: sll_t_particle_mesh_coupling_spline_cl_1d
     type(sll_t_spline_pp_1d) :: spline_pp
     
     ! Information about the particles
     sll_int32  :: no_particles 
     sll_int32  :: n_span 
     sll_real64 :: scaling 
     sll_int32  :: n_quad_points 
     sll_real64, allocatable :: spline_val(:) 
     sll_real64, allocatable :: spline_val_more(:) 
     ! For 1d quadrature 
     sll_int32 :: n_quad_points_line 
     sll_real64, allocatable :: quad_xw_line(:,:) 
     sll_real64, allocatable :: quad_xw(:,:) 
 
   contains
     procedure :: add_charge => add_charge_single_spline_cl_1d
     procedure :: add_particle_mass => add_particle_mass_spline_cl_1d
     procedure :: add_particle_mass_full => add_particle_mass_spline_cl_1d_full
     procedure :: evaluate => evaluate_field_single_spline_cl_1d
     procedure :: evaluate_multiple => evaluate_multiple_spline_cl_1d
     procedure :: add_current_update_v => add_current_update_v_spline_cl_1d
     procedure :: add_current => add_current_spline_cl_1d
     procedure :: add_current_evaluate => add_current_evaluate_spline_cl_1d
 
     procedure :: init => init_spline_cl_1d
     procedure :: free => free_spline_cl_1d
 
  end type sll_t_particle_mesh_coupling_spline_cl_1d
 
contains
 
 
  !---------------------------------------------------------------------------!
  subroutine add_charge_single_spline_cl_1d(self, position, marker_charge, rho_dofs)
    class( sll_t_particle_mesh_coupling_spline_cl_1d ), intent(inout)   :: self
    sll_real64,                               intent( in )    :: position(self%dim) 
    sll_real64,                               intent( in )    :: marker_charge 
    sll_real64,                               intent( inout ) :: rho_dofs(self%n_dofs) 
    !local variables
    sll_int32 :: i
    sll_int32 :: box
    sll_real64 :: xi
 
    call convert_x_to_xbox( self, position, xi, box )
    call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp, self%n_cells, self%spline_degree, xi, box, self%spline_val )
 
    do i = 1, self%n_span
       rho_dofs(box+i-1) = rho_dofs(box+i-1) +&
            (marker_charge * self%spline_val(i)* self%scaling)
    end do
 
 
  end subroutine add_charge_single_spline_cl_1d
 
 
  !---------------------------------------------------------------------------!
  subroutine add_particle_mass_spline_cl_1d(self, position, marker_charge, particle_mass)
    class( sll_t_particle_mesh_coupling_spline_cl_1d ), intent(inout)   :: self
    sll_real64,                               intent( in )    :: position(self%dim) 
    sll_real64,                               intent( in )    :: marker_charge 
    sll_real64,                               intent( inout ) :: particle_mass(:,:) 
 
    !local variables
    sll_int32 :: i, column
    sll_int32 :: box
    sll_real64 :: xi
 
    sll_assert( size(particle_mass,1) == self%n_span )
    sll_assert( size(particle_mass,2) == self%n_dofs )
 
    call convert_x_to_xbox( self, position, xi, box )
 
    call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp, self%n_cells, self%spline_degree, xi, box, self%spline_val )
 
    do i = 1, self%n_span
       do column = i, self%n_span
          particle_mass(column-i+1, box+i-1) = particle_mass( column-i+1, box+i-1) + &
               marker_charge * self%spline_val(i) * self%spline_val(column)
          ! Note: No scaling since Galerkin function
       end do
    end do
 
  end subroutine add_particle_mass_spline_cl_1d
 
 
  !----------------------------------------------------------------------!
  subroutine add_particle_mass_spline_cl_1d_full(self, position, marker_charge, particle_mass)
    class( sll_t_particle_mesh_coupling_spline_cl_1d ), intent(inout)   :: self
    sll_real64,                               intent( in )    :: position(self%dim) 
    sll_real64,                               intent( in )    :: marker_charge 
    sll_real64,                               intent( inout ) :: particle_mass(:,:) 
    !local variables
    sll_int32 :: i, column, ind
    sll_int32 :: box
    sll_real64 :: xi
 
    sll_assert( size(particle_mass,1) == 2*self%spline_degree+1 )
    sll_assert( size(particle_mass,2) == self%n_dofs )
 
    call convert_x_to_xbox( self, position, xi, box )
 
    call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp, self%n_cells, self%spline_degree, xi, box, self%spline_val )
 
    do i = 1, self%n_span
       ind=1+(self%n_span-i)
       do column = 1, self%n_span
          particle_mass(ind, box+i-1) = particle_mass( ind, box+i-1) + &
               marker_charge * self%spline_val(i) * self%spline_val(column)
          ind = ind+1
       end do
    end do
 
  end subroutine add_particle_mass_spline_cl_1d_full
 
 
  !---------------------------------------------------------------------------!
  subroutine evaluate_field_single_spline_cl_1d(self, position, field_dofs, field_value)
    class(sll_t_particle_mesh_coupling_spline_cl_1d), intent( inout ) :: self 
    sll_real64,                              intent( in )    :: position(self%dim) 
    sll_real64,                              intent( in )    :: field_dofs(self%n_dofs) 
    sll_real64,                              intent( out )   :: field_value 
    !local variables
    sll_int32 :: i
    sll_int32 :: box
    sll_real64 :: xi
 
    call convert_x_to_xbox( self, position, xi, box )
 
    call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp, self%n_cells, self%spline_degree, xi, box, self%spline_val )
 
 
    field_value = 0.0_f64
    do i = 1, self%n_span
       field_value = field_value + &
            field_dofs(box+i-1) *  &
            self%spline_val(i)
    end do
 
  end subroutine evaluate_field_single_spline_cl_1d
 
 
  !---------------------------------------------------------------------------!
  subroutine evaluate_multiple_spline_cl_1d(self, position, components, field_dofs, field_value)
    class(sll_t_particle_mesh_coupling_spline_cl_1d), intent( inout ) :: self 
    sll_real64,                              intent( in )    :: position(self%dim) 
    sll_int32,                               intent(in)      :: components(:) 
    sll_real64,                              intent( in )    :: field_dofs(:,:) 
    sll_real64,                              intent(out)     :: field_value(:) 
    !local variables
    sll_int32 :: i
    sll_int32 :: box
    sll_real64 :: xi
 
    sll_assert( size(field_dofs,1) == self%n_dofs )
    sll_assert( size(field_dofs,2) == size(field_value) )
 
    call convert_x_to_xbox( self, position, xi, box )
 
    call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp, self%n_cells, self%spline_degree, xi, box, self%spline_val )
 
 
    field_value = 0.0_f64
    do i = 1, self%n_span
 
       field_value = field_value + &
            field_dofs(box+i-1,components) *  &
            self%spline_val(i)
    end do
 
  end subroutine evaluate_multiple_spline_cl_1d
 
 
  !---------------------------------------------------------------------------!
  subroutine add_current_spline_cl_1d( self, position_old, position_new, marker_charge, j_dofs )
    class( sll_t_particle_mesh_coupling_spline_cl_1d ), intent(inout)   :: self
    sll_real64,                               intent( in )    :: position_old(self%dim) 
    sll_real64,                               intent( in )    :: position_new(self%dim) 
    sll_real64,                               intent( in )    :: marker_charge 
    sll_real64,                               intent( inout ) :: j_dofs(self%n_dofs) 
    !local variables
    sll_int32  :: index_old, index_new, ind
    sll_real64 :: r_old, r_new
 
    call convert_x_to_xbox( self, position_old, r_old, index_old )
    call convert_x_to_xbox( self, position_new, r_new, index_new )
 
      if (index_old == index_new) then
       if (r_old < r_new) then
          call update_j( self, r_old, r_new, index_old, marker_charge, &
               1.0_f64, j_dofs)
       else
          call update_j( self, r_new, r_old, index_old, marker_charge, &
               -1.0_f64, j_dofs)
       end if
    elseif (index_old < index_new) then
       call update_j ( self, r_old, 1.0_f64, index_old, marker_charge, &
            1.0_f64, j_dofs)
       call update_j ( self, 0.0_f64, r_new, index_new, marker_charge, &
            1.0_f64, j_dofs)
       do ind = index_old+1, index_new-1
          call update_j ( self, 0.0_f64, 1.0_f64, ind, marker_charge, &
               1.0_f64, j_dofs)
       end do
    else
       call update_j ( self, r_new, 1.0_f64, index_new, marker_charge, &
            -1.0_f64, j_dofs)
       call update_j ( self, 0.0_f64, r_old, index_old, marker_charge, &
            -1.0_f64, j_dofs)
       do ind = index_new+1, index_old-1
          call update_j ( self, 0.0_f64, 1.0_f64, ind, marker_charge, &
               -1.0_f64, j_dofs)
       end do
    end if
 
  end subroutine add_current_spline_cl_1d
 
    
  subroutine update_j(self, lower, upper, index, marker_charge,sign, j_dofs)
    class(sll_t_particle_mesh_coupling_spline_cl_1d), intent(inout) :: self
    sll_real64,                             intent(in)    :: lower
    sll_real64,                             intent(in)    :: upper
    sll_int32,                              intent(in)    :: index
    sll_real64,                             intent(in)    :: marker_charge
    sll_real64,                             intent(in)    :: sign
    sll_real64,                             intent(inout) :: j_dofs(:)
    !Local variables
    sll_int32  :: ind, i_grid, j
    sll_real64 :: c1, c2
 
    c1 =  0.5_f64*(upper-lower)
    c2 =  0.5_f64*(upper+lower)
 
    call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp, self%n_cells, self%spline_degree, c1*self%quad_xw(1,1)+c2, index, self%spline_val )
    self%spline_val = self%spline_val * (self%quad_xw(2,1)*c1)
    do j=2,self%n_quad_points
       call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp, self%n_cells, self%spline_degree, c1*self%quad_xw(1,j)+c2, index, self%spline_val_more )
       self%spline_val = self%spline_val + self%spline_val_more * (self%quad_xw(2,j)*c1)
    end do
    self%spline_val = self%spline_val * (sign*self%delta_x)
 
    ind = 1
    do i_grid = index, index + self%spline_degree
       j_dofs(i_grid) = j_dofs(i_grid) + &
            (marker_charge*self%spline_val(ind)* self%scaling)
       ind = ind + 1
    end do
 
  end subroutine update_j
  
  
  subroutine add_current_evaluate_spline_cl_1d(self, position_old, position_new, marker_charge, vbar, field_dofs, j_dofs, field)
    class( sll_t_particle_mesh_coupling_spline_cl_1d ), intent(inout)   :: self
    sll_real64,                               intent( in )    :: position_old(self%dim) 
    sll_real64,                               intent( in )    :: position_new(self%dim) 
    sll_real64,                               intent( in )    :: marker_charge 
    sll_real64,                               intent( in )    :: vbar 
    sll_real64,                               intent( in    ) :: field_dofs(self%n_dofs) 
    sll_real64,                               intent( inout ) :: j_dofs(self%n_dofs) 
    sll_real64,                               intent( out   ) :: field 
    !local variables
    sll_int32  :: index_old, index_new, ind
    sll_real64 :: r_old, r_new
 
    call convert_x_to_xbox( self, position_old, r_old, index_old )
    call convert_x_to_xbox( self, position_new, r_new, index_new )
 
    field = 0.0_f64
    
    if (index_old == index_new) then
       if (r_old < r_new) then
          call update_je( self, r_old, r_new, index_old, marker_charge, &
               1.0_f64, j_dofs, field_dofs, field)
       else
          call update_je( self, r_new, r_old, index_old, marker_charge, &
               -1.0_f64, j_dofs, field_dofs, field)
       end if
    elseif (index_old < index_new) then
       call update_je( self, r_old, 1.0_f64, index_old, marker_charge, &
            1.0_f64, j_dofs, field_dofs, field)
       call update_je( self, 0.0_f64, r_new, index_new, marker_charge, &
            1.0_f64, j_dofs, field_dofs, field)
       do ind = index_old+1, index_new-1
          call update_je( self, 0.0_f64, 1.0_f64, ind, marker_charge, &
               1.0_f64, j_dofs, field_dofs, field)
       end do
    else
       call update_je( self, r_new, 1.0_f64, index_new, marker_charge, &
            -1.0_f64, j_dofs, field_dofs, field)
       call update_je( self, 0.0_f64, r_old, index_old, marker_charge, &
            -1.0_f64, j_dofs, field_dofs, field)
       do ind = index_new+1, index_old-1
          call update_je( self, 0.0_f64, 1.0_f64, ind, marker_charge, &
               -1.0_f64, j_dofs, field_dofs, field)
       end do
    end if
    field = field/vbar
 
  end subroutine add_current_evaluate_spline_cl_1d
 
  
  subroutine update_je(self, lower, upper, index, marker_charge, sign, j_dofs, field_dofs, field)
    class(sll_t_particle_mesh_coupling_spline_cl_1d), intent(inout) :: self
    sll_real64,                             intent(in)    :: lower
    sll_real64,                             intent(in)    :: upper
    sll_int32,                              intent(in)    :: index
    sll_real64,                             intent(in)    :: marker_charge
    sll_real64,                             intent(in)    :: sign
    sll_real64,                             intent(inout) :: j_dofs(:)
    sll_real64,                             intent(in)    :: field_dofs(:)
    sll_real64,                             intent(inout)    :: field
    !Local variables
    sll_int32  :: ind, i_grid, j
    sll_real64 :: c1, c2
 
    c1 =  0.5_f64*(upper-lower)
    c2 =  0.5_f64*(upper+lower)
 
    call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp, self%n_cells, self%spline_degree, c1*self%quad_xw(1,1)+c2, index, self%spline_val )
    self%spline_val = self%spline_val * (self%quad_xw(2,1)*c1)
    do j=2,self%n_quad_points
       call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp, self%n_cells, self%spline_degree, c1*self%quad_xw(1,j)+c2, index, self%spline_val_more )
       self%spline_val = self%spline_val + self%spline_val_more * (self%quad_xw(2,j)*c1)
    end do
    self%spline_val = self%spline_val * self%delta_x * sign
 
    ind = 1
    do i_grid = index, index + self%spline_degree
       j_dofs(i_grid) = j_dofs(i_grid) + &
            (marker_charge*self%spline_val(ind)* self%scaling)
       field = field + self%spline_val(ind)*field_dofs(i_grid)
       ind = ind + 1
    end do
 
  end subroutine update_je
 
  
  !---------------------------------------------------------------------------!
  subroutine add_current_update_v_spline_cl_1d (self, position_old, position_new, marker_charge, qoverm, bfield_dofs, vi, j_dofs)
    class(sll_t_particle_mesh_coupling_spline_cl_1d), intent(inout) :: self
    sll_real64, intent(in)    :: position_old(self%dim) 
    sll_real64, intent(in)    :: position_new(self%dim) 
    sll_real64, intent(in)    :: marker_charge 
    sll_real64, intent(in)    :: qoverm 
    sll_real64, intent(in)    :: bfield_dofs(:) 
    sll_real64, intent(inout) :: vi(:) 
    sll_real64, intent(inout) :: j_dofs(:) 
    ! local variables
    sll_int32  :: index_old, index_new, ind
    sll_real64 :: r_old, r_new
 
    call convert_x_to_xbox( self, position_old, r_old, index_old )
    call convert_x_to_xbox( self, position_new, r_new, index_new )
 
    if (index_old == index_new) then
       if (r_old < r_new) then
          call update_jv( self, r_old, r_new, index_old, marker_charge, &
               qoverm, 1.0_f64, vi(2), j_dofs, bfield_dofs)
       else
          call update_jv( self, r_new, r_old, index_old, marker_charge, qoverm, &
               -1.0_f64, vi(2), j_dofs, bfield_dofs)
       end if
    elseif (index_old < index_new) then
       call update_jv ( self, r_old, 1.0_f64, index_old, marker_charge, &
            qoverm, 1.0_f64, vi(2), j_dofs, bfield_dofs)
       call update_jv ( self, 0.0_f64, r_new, index_new, marker_charge, &
            qoverm, 1.0_f64, vi(2), j_dofs, bfield_dofs)
       do ind = index_old+1, index_new-1
          call update_jv ( self, 0.0_f64, 1.0_f64, ind, marker_charge, &
               qoverm, 1.0_f64, vi(2), j_dofs, bfield_dofs)
       end do
    else
       call update_jv ( self, r_new, 1.0_f64, index_new, marker_charge, qoverm, &
            -1.0_f64, vi(2), j_dofs, bfield_dofs)
       call update_jv ( self, 0.0_f64, r_old, index_old, marker_charge, qoverm, &
            -1.0_f64, vi(2), j_dofs, bfield_dofs)
       do ind = index_new+1, index_old-1
          call update_jv ( self, 0.0_f64, 1.0_f64, ind, marker_charge, qoverm, &
               -1.0_f64, vi(2), j_dofs, bfield_dofs)
       end do
    end if
 
  end subroutine add_current_update_v_spline_cl_1d
 
  
  subroutine update_jv(self, lower, upper, index, marker_charge, qoverm, sign, vi, j_dofs, bfield_dofs)
    class(sll_t_particle_mesh_coupling_spline_cl_1d), intent(inout) :: self
    sll_real64,                             intent(in)    :: lower
    sll_real64,                             intent(in)    :: upper
    sll_int32,                              intent(in)    :: index
    sll_real64,                             intent(in)    :: marker_charge
    sll_real64,                             intent(in)    :: qoverm
    sll_real64,                             intent(in)    :: sign
    sll_real64,                             intent(inout) :: vi
    sll_real64,                             intent(in)    :: bfield_dofs(:)
    sll_real64,                             intent(inout) :: j_dofs(:)
    !Local variables
    sll_int32  :: ind, i_grid, j
    sll_real64 :: c1, c2
 
    c1 =  0.5_f64*(upper-lower)
    c2 =  0.5_f64*(upper+lower)
 
    call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp, self%n_cells, self%spline_degree, c1*self%quad_xw(1,1)+c2, index, self%spline_val )
    self%spline_val = self%spline_val * (self%quad_xw(2,1)*c1)
    do j=2,self%n_quad_points
       call sll_s_uniform_bsplines_eval_basis_clamped( self%spline_pp, self%n_cells, self%spline_degree, c1*self%quad_xw(1,j)+c2, index, self%spline_val_more )
       self%spline_val = self%spline_val + self%spline_val_more * (self%quad_xw(2,j)*c1)
    end do
    self%spline_val = self%spline_val * (sign*self%delta_x)
 
    ind = 1
    do i_grid = index, index + self%spline_degree
       j_dofs(i_grid) = j_dofs(i_grid) + &
            (marker_charge*self%spline_val(ind)* self%scaling)
       vi = vi - qoverm* self%spline_val(ind)*bfield_dofs(i_grid)
       ind = ind + 1
    end do
 
  end subroutine update_jv
 
 
  
  subroutine convert_x_to_xbox( self, position, xi, box )
    class( sll_t_particle_mesh_coupling_spline_cl_1d ), intent(inout)   :: self
    sll_real64,                               intent( in )    :: position(self%dim) 
    sll_real64,                               intent( out )    :: xi 
    sll_int32,                                intent( out )    :: box 
 
    xi = (position(1) - self%domain(1)) /self%delta_x
    box = floor( xi ) + 1
    xi = xi - real(box-1, f64)
 
    if( box == self%n_cells + 1 ) then
       if( xi == 0._f64)then 
          xi = 1._f64
          box = self%n_cells
       else
          print*, 'box,x, xbox', box, position(1), xi
          sll_error('convert_x_to_xbox', 'box1 value to high' ) 
       end if
    else if( box > self%n_cells + 1 ) then
       print*, 'box,x, xbox', box, position(1), xi
       sll_error('convert_x_to_xbox', 'box1 value to high' ) 
    else if( box <= 0 ) then
       print*, 'box,x, xbox', box, position(1), xi
       sll_error('convert_x_to_xbox', 'box1 value to low' ) 
    end if
 
  end subroutine convert_x_to_xbox
 
  subroutine sll_s_uniform_bsplines_eval_basis_clamped( spline, n_cells, degree, xi, box, spline_val )
    type(sll_t_spline_pp_1d), intent( in    ) :: spline
    sll_int32,                intent( in    ) :: n_cells
    sll_int32,                intent( in    ) :: degree
    sll_real64,               intent( in    ) :: xi
    sll_int32,                intent( in    ) :: box
    sll_real64,               intent(   out ) :: spline_val(:)
    !local variables
    sll_int32 :: i
 
    if(box >= 1 .and. box <= degree-1)then
       do i=1, degree+1
          spline_val(i) = sll_f_spline_pp_horner_1d( degree, spline%poly_coeffs_boundary_left(:,:,box), xi, i)
       end do
    else if (box >= degree .and. box <= n_cells-degree+1)then
       call sll_s_uniform_bsplines_eval_basis( degree, xi, spline_val )
    else if(box >= n_cells-degree+2 .and. box <= n_cells)then
       do i=1, degree+1
          spline_val(i) = sll_f_spline_pp_horner_1d( degree, spline%poly_coeffs_boundary_right(:,:,box-1-n_cells+degree), xi, i)
       end do
    else
       print*, 'box, xbox', box, xi
       sll_error( "uniform_bsplines_eval_basis_clamped", "box out of range" )
    end if
 
  end subroutine sll_s_uniform_bsplines_eval_basis_clamped
 
 
  subroutine init_spline_cl_1d( self, domain, n_cells, no_particles, spline_degree, smoothing_type, boundary )
    class( sll_t_particle_mesh_coupling_spline_cl_1d), intent(out)  :: self
    sll_real64,                              intent(in) :: domain(2) 
    sll_int32,                               intent(in) :: n_cells 
    sll_int32,                               intent(in) :: no_particles 
    sll_int32,                               intent(in) :: spline_degree 
    sll_int32,                               intent(in) :: smoothing_type 
    sll_int32,                               intent(in) :: boundary 
    !local variables
    sll_int32 :: ierr
 
 
    self%dim = 1
 
    ! Store grid information
    self%domain = domain
    self%n_cells = n_cells
    self%n_dofs = n_cells+spline_degree
    self%delta_x = (self%domain(2)-self%domain(1))/real(n_cells, f64)
 
    ! Store basis function information
    self%no_particles = no_particles
 
    ! Initialize information on the spline
    self%spline_degree = spline_degree
    self%n_span = spline_degree + 1
 
    ! Initialize information on smoothing type
    if (smoothing_type == sll_p_collocation) then
       self%scaling = 1.0_f64/self%delta_x
    elseif (smoothing_type == sll_p_galerkin) then
       self%scaling = 1.0_f64
    else
       print*, 'Smoothing Type ', smoothing_type, ' not implemented for kernel_smoother_spline_1d.'
    end if
 
    self%n_quad_points = (self%spline_degree+2)/2
 
    ALLOCATE( self%spline_val(self%n_span), stat = ierr)
    sll_assert( ierr == 0 )
    ALLOCATE( self%spline_val_more(self%n_span), stat = ierr )
    sll_assert( ierr == 0 )
    ALLOCATE( self%quad_xw(2,self%n_quad_points), stat = ierr )
    sll_assert( ierr == 0 )
    
    ! normalized Gauss Legendre points and weights
    self%quad_xw = sll_f_gauss_legendre_points_and_weights(self%n_quad_points)
 
    call sll_s_spline_pp_init_1d( self%spline_pp, spline_degree, n_cells, boundary)
 
    ! for line integral
    self%n_quad_points_line = (self%spline_degree+2)/2
 
    allocate( self%quad_xw_line(2,self%n_quad_points_line) )
    ! normalized Gauss Legendre points and weights
    self%quad_xw_line = sll_f_gauss_legendre_points_and_weights(self%n_quad_points_line)
 
    self%quad_xw_line(1,:) = 0.5_f64*(self%quad_xw_line(1,:)+1.0_f64)
    self%quad_xw_line(2,:) = 0.5_f64*self%quad_xw_line(2,:)
 
 
  end subroutine init_spline_cl_1d
 
 
  subroutine free_spline_cl_1d(self)
    class(sll_t_particle_mesh_coupling_spline_cl_1d), intent( inout ) :: self 
 
    deallocate(self%spline_val)
    call sll_s_spline_pp_free_1d(self%spline_pp)
 
  end subroutine free_spline_cl_1d
 
 
  !-------------------------------------------------------------------------------------------
  
  subroutine sll_s_new_particle_mesh_coupling_spline_cl_1d_ptr(smoother, domain, n_cells, no_particles, spline_degree, smoothing_type, boundary)
    class( sll_c_particle_mesh_coupling_1d), pointer, intent(out) :: smoother
    sll_int32,                              intent(in)  :: n_cells 
    sll_real64,                             intent(in)  :: domain(2) 
    sll_int32,                              intent(in)  :: no_particles 
    sll_int32,                              intent(in)  :: spline_degree 
    sll_int32,                              intent(in)  :: smoothing_type 
    sll_int32,                               intent(in) :: boundary 
    !local variables
    sll_int32 :: ierr
 
 
    sll_allocate( sll_t_particle_mesh_coupling_spline_cl_1d :: smoother , ierr)
    sll_assert( ierr == 0)
 
    select type( smoother )
    type is ( sll_t_particle_mesh_coupling_spline_cl_1d )
       call smoother%init( domain, n_cells, no_particles, spline_degree, smoothing_type, boundary )
    end select
 
  end subroutine sll_s_new_particle_mesh_coupling_spline_cl_1d_ptr
 
 
  !----------------------------------------------------------------------------------
  
  subroutine sll_s_new_particle_mesh_coupling_spline_cl_1d(smoother, domain, n_cells, no_particles, spline_degree, smoothing_type, boundary)
    class( sll_c_particle_mesh_coupling_1d), allocatable, intent(out):: smoother
    sll_int32,                                  intent(in) :: n_cells 
    sll_real64,                                 intent(in) :: domain(2) 
    sll_int32,                                  intent(in) :: no_particles 
    sll_int32,                                  intent(in) :: spline_degree 
    sll_int32,                                  intent(in) :: smoothing_type 
    sll_int32,                                  intent(in) :: boundary 
    !local variables
    sll_int32 :: ierr
 
 
    allocate( sll_t_particle_mesh_coupling_spline_cl_1d :: smoother , stat=ierr)
    sll_assert( ierr == 0)
 
    select type( smoother )
    type is ( sll_t_particle_mesh_coupling_spline_cl_1d )
       call smoother%init( domain, n_cells, no_particles, spline_degree, smoothing_type, boundary )
    end select
 
  end subroutine sll_s_new_particle_mesh_coupling_spline_cl_1d
 
 
end module sll_m_particle_mesh_coupling_spline_cl_1d