sll_m_time_propagator_pic_vm_1d2v_boris.F90

Go to the documentation of this file.

 
module sll_m_time_propagator_pic_vm_1d2v_boris
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_assert.h"
#include "sll_memory.h"
#include "sll_working_precision.h"
 
  use sll_m_collective, only: &
    sll_o_collective_allreduce, &
    sll_v_world_collective
 
  use sll_m_time_propagator_base, only: &
    sll_c_time_propagator_base
 
  use sll_m_particle_mesh_coupling_base_1d, only: &
    sll_c_particle_mesh_coupling_1d
 
  use sll_m_maxwell_1d_base, only: &
    sll_c_maxwell_1d_base
 
  use sll_m_particle_group_base, only: &
    sll_t_particle_array
 
  use mpi, only: &
    mpi_sum
 
  implicit none
 
  public :: &
    sll_t_time_propagator_pic_vm_1d2v_boris
 
  private
 
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
  type, extends(sll_c_time_propagator_base) :: sll_t_time_propagator_pic_vm_1d2v_boris
     class(sll_c_maxwell_1d_base), pointer :: maxwell_solver
     class(sll_c_particle_mesh_coupling_1d), pointer :: kernel_smoother_0
     class(sll_c_particle_mesh_coupling_1d), pointer :: kernel_smoother_1
     class(sll_t_particle_array), pointer  :: particle_group
         sll_int32 :: spline_degree 
     sll_real64 :: lx 
     sll_real64 :: x_min 
     sll_real64 :: delta_x 
 
     sll_real64 :: cell_integrals_0(4) 
     sll_real64 :: cell_integrals_1(3) 
 
 
     sll_real64, pointer     :: efield_dofs(:,:) 
     sll_real64, pointer     :: efield_dofs_mid(:,:) 
     sll_real64, pointer     :: bfield_dofs(:)   
     sll_real64, allocatable :: bfield_dofs_mid(:) 
     sll_real64, allocatable :: j_dofs(:,:)      
     sll_real64, allocatable :: j_dofs_local(:,:)
 
   contains
     procedure :: lie_splitting => lie_boris
     procedure :: lie_splitting_back => lie_boris
     procedure :: strang_splitting => operator_boris
 
     procedure :: init => initialize_pic_vm_1d2v_boris
     procedure :: free => delete_pic_vm_1d2v_boris
 
     procedure :: staggering => staggering_pic_vm_1d2v_boris
 
  end type sll_t_time_propagator_pic_vm_1d2v_boris
 
contains
 
  subroutine lie_boris(self, dt, number_steps)
    class(sll_t_time_propagator_pic_vm_1d2v_boris), intent(inout) :: self
    sll_real64,      intent(in)    :: dt   
    sll_int32,                                      intent(in)    :: number_steps 
 
 
    print*, 'Lie splitting not implemented in sll_m_time_propagator_pic_vm_1d2v_boris.'
 
  end subroutine lie_boris
 
  subroutine operator_boris(self, dt, number_steps)
    class(sll_t_time_propagator_pic_vm_1d2v_boris), intent(inout) :: self
    sll_real64,                                     intent(in)    :: dt   
    sll_int32,                                      intent(in)    :: number_steps 
 
    sll_int32 :: i_step
 
    do i_step = 1, number_steps
       ! (1) Compute B_{n+1/2} from B_n
       self%bfield_dofs_mid = self%bfield_dofs
       call self%maxwell_solver%compute_B_from_E( &
            dt, self%efield_dofs_mid(:,2), self%bfield_dofs)
       self%bfield_dofs_mid = (self%bfield_dofs_mid + self%bfield_dofs)*0.5_f64
       
       ! (2) Propagate v: v_{n-1/2} -> v_{n+1/2}
       ! (2a) Half time step with E-part
       call push_v_epart( self, dt*0.5_f64 )
       ! (2b) Full time step with B-part
       call push_v_bpart( self, dt )
       ! (2c) Half time step with E-part
       call push_v_epart( self, dt*0.5_f64 )
       
       ! (3) Propagate x: x_n -> x_{n+1}. Includes also accumulation of j_x, j_y
       call push_x_accumulate_j ( self, dt )
       
       ! (4) Compute E_{n+1}       
       self%efield_dofs = self%efield_dofs_mid
       ! Ex part:
       call self%maxwell_solver%compute_E_from_j(self%j_dofs(:,1), 1, self%efield_dofs_mid(:,1))
       ! TODO: Combine the two steps for efficiency
       ! Ey part:    
       call self%maxwell_solver%compute_E_from_B(&
            dt, self%bfield_dofs, self%efield_dofs_mid(:,2))
       call self%maxwell_solver%compute_E_from_j(self%j_dofs(:,2), 2, self%efield_dofs_mid(:,2))
       
    end do
 
  end subroutine operator_boris
 
  subroutine push_v_epart (self, dt)
    class(sll_t_time_propagator_pic_vm_1d2v_boris), intent(inout) :: self
    sll_real64,                                     intent(in)    :: dt   
 
    !local variables
    sll_int32 :: i_part, i_sp
    sll_real64 :: v_new(3), xi(3)
    sll_real64 :: efield(2)
    sll_real64 :: qm
 
 
    do i_sp =1, self%particle_group%n_species
       qm = self%particle_group%group(i_sp)%species%q_over_m();
       ! V_new = V_old + dt * E
       do i_part=1,self%particle_group%group(i_sp)%n_particles
          ! Evaluate efields at particle position
          xi = self%particle_group%group(i_sp)%get_x(i_part)
          call self%kernel_smoother_1%evaluate &
               (xi(1), self%efield_dofs_mid(:,1), efield(1))
          call self%kernel_smoother_0%evaluate &
               (xi(1), self%efield_dofs_mid(:,2), efield(2))
          v_new = self%particle_group%group(i_sp)%get_v(i_part)
          v_new(1:2) = v_new(1:2) + dt* qm * efield
          call self%particle_group%group(i_sp)%set_v(i_part, v_new)
       end do
    end do
    
 
  end subroutine push_v_epart
 
   subroutine push_v_bpart (self, dt)
    class(sll_t_time_propagator_pic_vm_1d2v_boris), intent(inout) :: self
    sll_real64,                                     intent(in)    :: dt   
 
    !local variables
    sll_int32 :: i_part, i_sp
    sll_real64 :: vi(3), v_new(3), xi(3)
    sll_real64 :: bfield, m11, m12
    sll_real64 :: qmdt
 
    do i_sp =1, self%particle_group%n_species
       qmdt = self%particle_group%group(i_sp)%species%q_over_m()*0.5_f64*dt;
       
       do i_part=1,self%particle_group%group(i_sp)%n_particles
          vi= self%particle_group%group(i_sp)%get_v(i_part)
          xi = self%particle_group%group(i_sp)%get_x(i_part)
          call self%kernel_smoother_1%evaluate &
               (xi(1), self%bfield_dofs_mid, bfield)
          
          bfield = qmdt*bfield
          m11 = 1.0_f64/(1.0_f64 + bfield**2) 
          m12 = m11*bfield*2.0_f64
          m11 = m11*(1-bfield**2)
          
          v_new(1) = m11 * vi(1) + m12 * vi(2)
          v_new(2) = - m12 * vi(1) + m11 * vi(2)
          v_new(3) = 0.0_f64
          
          call self%particle_group%group(i_sp)%set_v(i_part, v_new)
 
       end do
    end do
 
  end subroutine push_v_bpart
 
  subroutine push_x_accumulate_j (self, dt)
    class(sll_t_time_propagator_pic_vm_1d2v_boris), intent(inout) :: self
    sll_real64,                                     intent(in)    :: dt   
 
    !local variables
    sll_int32 :: i_part, i_sp
    sll_real64 :: x_new(3), vi(3), wi(1), x_old(3)
    sll_int32  :: n_cells
    sll_real64 :: qoverm
 
 
    n_cells = self%kernel_smoother_0%n_dofs
 
 
    self%j_dofs_local = 0.0_f64
 
    ! For each particle compute the index of the first DoF on the grid it contributes to and its position (normalized to cell size one). Note: j_dofs(_local) does not hold the values for j itself but for the integrated j.
    ! Then update particle position:  X_new = X_old + dt * V
 
    
    do i_sp =1, self%particle_group%n_species
       qoverm = self%particle_group%group(i_sp)%species%q_over_m();
       do i_part=1,self%particle_group%group(i_sp)%n_particles  
          ! Read out particle position and velocity
          x_old = self%particle_group%group(i_sp)%get_x(i_part)
          vi = self%particle_group%group(i_sp)%get_v(i_part)
          
          ! Then update particle position:  X_new = X_old + dt * V
          x_new = x_old + dt * vi
          
          ! Get charge for accumulation of j
          wi = self%particle_group%group(i_sp)%get_charge(i_part)
          
          ! TODO: Check here also second posibility with sum of two accumulations
!!$          ! Accumulate jx
!!$          call self%kernel_smoother_1%add_charge( [(x_old(1)+x_new(1))*0.5_f64], wi(1)*vi(1)*dt, &
!!$               self%j_dofs_local(:,1))
!!$          ! Accumulate jy
!!$          call self%kernel_smoother_0%add_charge( [(x_old(1)+x_new(1))*0.5_f64], wi(1)*vi(2)*dt, &
!!$               self%j_dofs_local(:,2))
          call self%kernel_smoother_1%add_current( x_old(1), x_new(1), wi(1), self%j_dofs_local(:,1) )
          if ( abs(vi(1)) > 1e-15_f64 ) then
             call self%kernel_smoother_0%add_current( x_old(1), x_new(1), wi(1)*vi(2)/vi(1), self%j_dofs_local(:,2) )
          else
             call self%kernel_smoother_0%add_charge( [(x_old(1)+x_new(1))*0.5_f64], wi(1)*vi(2)*dt , &
                  self%j_dofs_local(:,2))
          end if
          
          x_new(1) = modulo(x_new(1), self%Lx)
          call self%particle_group%group(i_sp)%set_x(i_part, x_new)
          
       end do
    end do
 
    self%j_dofs = 0.0_f64
    ! MPI to sum up contributions from each processor
    call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local(:,1), &
         n_cells, mpi_sum, self%j_dofs(:,1))
    call sll_o_collective_allreduce( sll_v_world_collective, self%j_dofs_local(:,2), &
         n_cells, mpi_sum, self%j_dofs(:,2))
    
 
  end subroutine push_x_accumulate_j
 
 !---------------------------------------------------------------------------!
  subroutine initialize_pic_vm_1d2v_boris(&
       self, &
       maxwell_solver, &
       kernel_smoother_0, &
       kernel_smoother_1, &
       particle_group, &
       efield_dofs, &
       bfield_dofs, &
       x_min, &
       Lx) 
    class(sll_t_time_propagator_pic_vm_1d2v_boris), intent(out) :: self
    class(sll_c_maxwell_1d_base), target,          intent(in)  :: maxwell_solver
    class(sll_c_particle_mesh_coupling_1d), target,          intent(in)  :: kernel_smoother_0
    class(sll_c_particle_mesh_coupling_1d), target,          intent(in)  :: kernel_smoother_1
    class(sll_t_particle_array), target,      intent(in)  :: particle_group
    sll_real64, target,                            intent(in)  :: efield_dofs(:,:) 
    sll_real64, target,                            intent(in)  :: bfield_dofs(:) 
    sll_real64,                                     intent(in)  :: x_min 
    sll_real64,                                     intent(in)  :: lx 
 
    !local variables
    sll_int32 :: ierr
 
    self%maxwell_solver => maxwell_solver
    self%kernel_smoother_0 => kernel_smoother_0
    self%kernel_smoother_1 => kernel_smoother_1
    self%particle_group => particle_group
    self%efield_dofs => efield_dofs
    self%bfield_dofs => bfield_dofs
 
#ifndef __PGI
    ! Check that n_dofs is the same for both kernel smoothers.
    sll_assert( self%kernel_smoother_0%n_dofs == self%kernel_smoother_1%n_dofs )
#endif
    sll_allocate(self%j_dofs(self%kernel_smoother_0%n_dofs,2), ierr)
    sll_allocate(self%j_dofs_local(self%kernel_smoother_0%n_dofs,2), ierr)
    sll_allocate(self%bfield_dofs_mid(self%kernel_smoother_1%n_dofs), ierr)
    sll_allocate(self%efield_dofs_mid(self%kernel_smoother_1%n_dofs,2), ierr)
 
    self%spline_degree = self%kernel_smoother_0%spline_degree
    self%x_min = x_min
    self%Lx = lx
    self%delta_x = self%Lx/real(self%kernel_smoother_1%n_dofs, f64)
    
    self%cell_integrals_1 = [0.5_f64, 2.0_f64, 0.5_f64]
    self%cell_integrals_1 = self%cell_integrals_1 / 3.0_f64
 
    self%cell_integrals_0 = [1.0_f64,11.0_f64,11.0_f64,1.0_f64]
    self%cell_integrals_0 = self%cell_integrals_0 / 24.0_f64
 
 
  end subroutine initialize_pic_vm_1d2v_boris
 
  !---------------------------------------------------------------------------!
  subroutine delete_pic_vm_1d2v_boris(self)
    class(sll_t_time_propagator_pic_vm_1d2v_boris), intent( inout ) :: self
 
    deallocate(self%j_dofs)
    deallocate(self%j_dofs_local)
    deallocate(self%bfield_dofs_mid)
    deallocate(self%efield_dofs_mid)
    self%maxwell_solver => null()
    self%kernel_smoother_0 => null()
    self%kernel_smoother_1 => null()
    self%particle_group => null()
    self%efield_dofs => null()
    self%bfield_dofs => null()
 
  end subroutine delete_pic_vm_1d2v_boris
 
 
  !---------------------------------------------------------------------------!
  subroutine staggering_pic_vm_1d2v_boris(self, dt)
    class(sll_t_time_propagator_pic_vm_1d2v_boris), intent( inout ) :: self
    sll_real64,                                     intent(in)    :: dt   
 
    call push_x_accumulate_j (self, dt*0.5_f64)
 
    ! (4) Compute E_{n+1}       
    self%efield_dofs_mid = self%efield_dofs
  !  self%j_dofs = dt*0.5_f64*self%j_dofs
    ! Ex part:
    call self%maxwell_solver%compute_E_from_j(self%j_dofs(:,1), 1, self%efield_dofs_mid(:,1))
    ! TODO: Combine the two steps for efficiency
    ! Ey part:    
    call self%maxwell_solver%compute_E_from_B(&
         dt*0.5_f64, self%bfield_dofs, self%efield_dofs_mid(:,2))
    call self%maxwell_solver%compute_E_from_j(self%j_dofs(:,2), 2, self%efield_dofs_mid(:,2))
 
  end subroutine staggering_pic_vm_1d2v_boris
 
 
end module sll_m_time_propagator_pic_vm_1d2v_boris