Lorentz_integrate_particles Subroutine

private subroutine Lorentz_integrate_particles(end_time)

Uses

  • proc~~lorentz_integrate_particles~~UsesGraph proc~lorentz_integrate_particles Lorentz_integrate_particles module~mod_geometry mod_geometry proc~lorentz_integrate_particles->module~mod_geometry module~mod_global_parameters mod_global_parameters proc~lorentz_integrate_particles->module~mod_global_parameters module~mod_usr_methods mod_usr_methods proc~lorentz_integrate_particles->module~mod_usr_methods module~mod_comm_lib mod_comm_lib module~mod_geometry->module~mod_comm_lib module~mod_basic_types mod_basic_types module~mod_global_parameters->module~mod_basic_types module~mod_connectivity mod_connectivity module~mod_global_parameters->module~mod_connectivity module~mod_constants mod_constants module~mod_global_parameters->module~mod_constants module~mod_physicaldata mod_physicaldata module~mod_global_parameters->module~mod_physicaldata module~mod_variables mod_variables module~mod_global_parameters->module~mod_variables mpi mpi module~mod_global_parameters->mpi module~mod_variables->module~mod_basic_types

Relativistic particle integrator

Arguments

Type IntentOptional Attributes Name
double precision, intent(in) :: end_time

Calls

proc~~lorentz_integrate_particles~~CallsGraph proc~lorentz_integrate_particles Lorentz_integrate_particles proc~fix_phi_crossing fix_phi_crossing proc~lorentz_integrate_particles->proc~fix_phi_crossing proc~get_bfield get_bfield proc~lorentz_integrate_particles->proc~get_bfield proc~get_efield get_efield proc~lorentz_integrate_particles->proc~get_efield proc~get_lfac get_lfac proc~lorentz_integrate_particles->proc~get_lfac proc~lorentz_get_particle_dt Lorentz_get_particle_dt proc~lorentz_integrate_particles->proc~lorentz_get_particle_dt proc~lorentz_kick Lorentz_kick proc~lorentz_integrate_particles->proc~lorentz_kick proc~lorentz_update_payload Lorentz_update_payload proc~lorentz_integrate_particles->proc~lorentz_update_payload proc~partcoord_from_cartesian partcoord_from_cartesian proc~lorentz_integrate_particles->proc~partcoord_from_cartesian proc~partcoord_to_cartesian partcoord_to_cartesian proc~lorentz_integrate_particles->proc~partcoord_to_cartesian proc~partvec_from_cartesian partvec_from_cartesian proc~lorentz_integrate_particles->proc~partvec_from_cartesian proc~partvec_to_cartesian partvec_to_cartesian proc~lorentz_integrate_particles->proc~partvec_to_cartesian proc~point_in_igrid_ghostc point_in_igrid_ghostc proc~fix_phi_crossing->proc~point_in_igrid_ghostc proc~get_vec get_vec proc~get_bfield->proc~get_vec proc~interpolate_var interpolate_var proc~get_bfield->proc~interpolate_var proc~get_efield->proc~get_bfield proc~get_efield->proc~get_vec proc~get_efield->proc~interpolate_var proc~lorentz_get_particle_dt->proc~get_lfac proc~lorentz_get_particle_dt->proc~get_vec proc~limit_dt_endtime limit_dt_endtime proc~lorentz_get_particle_dt->proc~limit_dt_endtime proc~lorentz_kick->proc~get_lfac proc~cross cross proc~lorentz_kick->proc~cross proc~lorentz_update_payload->proc~get_bfield proc~lorentz_update_payload->proc~get_efield proc~lorentz_update_payload->proc~get_lfac proc~lorentz_update_payload->proc~cross proc~mpistop mpistop proc~cross->proc~mpistop proc~get_vec->proc~interpolate_var proc~interpolate_var->proc~mpistop mpi_abort mpi_abort proc~mpistop->mpi_abort

Source Code

  subroutine Lorentz_integrate_particles(end_time)
    use mod_global_parameters
    use mod_geometry
    use mod_usr_methods, only: usr_create_particles, usr_update_payload
    double precision, intent(in)      :: end_time
    double precision                  :: defpayload(ndefpayload)
    double precision                  :: usrpayload(nusrpayload)
    integer                           :: ipart, iipart, igrid
    double precision                  :: lfac, q, m, dt_p
    double precision                  :: xp(ndir), xpm(ndir), xpc(ndir),&
        xpcm(ndir)
    double precision                  :: up(ndir), upc(ndir), tp
    double precision, dimension(ndir) :: b, e, bc, ec, vfluid, current
    double precision                  :: rho, rhoold, td

    do iipart=1,nparticles_active_on_mype
      ipart = particles_active_on_mype(iipart);
      q     = particle(ipart)%self%q
      m     = particle(ipart)%self%m
      igrid = particle(ipart)%igrid
      xp    = particle(ipart)%self%x
      up    = particle(ipart)%self%u
      tp    = particle(ipart)%self%time   

      dt_p = Lorentz_get_particle_dt(particle(ipart), end_time)
      particle(ipart)%self%dt = dt_p

      ! Push particle over first half time step
      ! all pushes done in Cartesian coordinates
      call partcoord_to_cartesian(xp,xpc)
      call partvec_to_cartesian(xp,up,upc)
      call get_lfac(upc,lfac)
      xpcm = xpc + 0.5d0 * dt_p * upc/lfac
      call partcoord_from_cartesian(xpm,xpcm)
      ! Fix xp if the 0,2*pi boundary was crossed in cylindrical/spherical coords
      call fix_phi_crossing(xpm,igrid)
      
      ! Get E, B at n+1/2 position
      call get_bfield(igrid, xpm, tp+dt_p/2.d0, b)
      call get_efield(igrid, xpm, tp+dt_p/2.d0, e)

!      call get_vec(vp, igrid, xpm, tp+dt_p/2.d0, vfluid)
!      call get_vec(jp, igrid, xpm, tp+dt_p/2.d0, current)
!      select case (coordinate)
!      case (Cartesian,Cartesian_stretched,spherical)
!        e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
!        e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
!        e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)
!      case (cylindrical)
!        e(r_) = -vfluid(phi_)*b(z_)+vfluid(z_)*b(phi_) + particles_eta*current(r_)
!        e(phi_) = vfluid(r_)*b(z_)-vfluid(z_)*b(r_) + particles_eta*current(phi_)
!        e(z_) = -vfluid(r_)*b(phi_)+vfluid(phi_)*b(r_) + particles_eta*current(z_)
!      end select
!      if (particles_etah > zero) then
!        call interpolate_var(igrid,ixG^LL,ixM^LL,ps(igrid)%w(ixG^T,1),ps(igrid)%x,xpm,rho)
!        if (time_advance) then
!          td = (tp+dt_p/2.d0 - global_time) / dt
!          call interpolate_var(igrid,ixG^LL,ixM^LL,pso(igrid)%w(ixG^T,1),ps(igrid)%x,xpm,rhoold)
!          rho = rhoold * (1.d0-td) + rho * td
!        end if
!        select case (coordinate)
!        case (Cartesian,Cartesian_stretched,spherical)
!          e(1) = e(1) + particles_etah/rho * (current(2)*b(3) - current(3)*b(2))
!          e(2) = e(2) + particles_etah/rho * (-current(1)*b(3) + current(3)*b(1))
!          e(3) = e(3) + particles_etah/rho * (current(1)*b(2) - current(2)*b(1))
!        case (cylindrical)
!          e(r_) = e(r_) + particles_etah/rho * (current(phi_)*b(z_) - current(z_)*b(phi_))
!          e(phi_) = e(phi_) + particles_etah/rho * (-current(r_)*b(z_) + current(z_)*b(r_))
!          e(z_) = e(z_) + particles_etah/rho * (current(r_)*b(phi_) - current(phi_)*b(r_))
!        end select
!      end if

      ! Convert fields to Cartesian frame
      call partvec_to_cartesian(xpm,b,bc)
      call partvec_to_cartesian(xpm,e,ec)

      ! 'Kick' particle (update velocity) based on the chosen integrator
      call Lorentz_kick(upc,ec,bc,q,m,dt_p)

      ! Push particle over second half time step
      ! all pushes done in Cartesian coordinates
      call get_lfac(upc,lfac)
      xpc = xpcm + 0.5d0 * dt_p * upc/lfac
      call partcoord_from_cartesian(xp,xpc)
      ! Fix xp if the 0,2*pi boundary was crossed in cylindrical/spherical coords
      call fix_phi_crossing(xp,igrid)
      call partvec_from_cartesian(xp,up,upc)

      ! Store updated x,u
      particle(ipart)%self%x = xp
      particle(ipart)%self%u = up

      ! Time update
      tp = tp + dt_p
      particle(ipart)%self%time = tp

      ! Update payload
      call Lorentz_update_payload(igrid,xp,up,q,m,defpayload,ndefpayload,tp)
      particle(ipart)%payload(1:ndefpayload) = defpayload
      if (associated(usr_update_payload)) then
        call usr_update_payload(igrid,xp,up,q,m,usrpayload,nusrpayload,tp)
        particle(ipart)%payload(ndefpayload+1:npayload) = usrpayload
      end if

    end do ! ipart loop

  end subroutine Lorentz_integrate_particles