mod_particle_gca.t

Go to the documentation of this file.

!> Particle mover with Newtonian/relativistic Guiding Center Approximation (GCA)
!> By Jannis Teunissen, Bart Ripperda, Oliver Porth, and Fabio Bacchini (2016-2020)
module mod_particle_gca
  use mod_particle_base
 
  private
 
  !> Variable index for gradient B, with relativistic correction 1/kappa
  !> where kappa = 1/sqrt(1 - E_perp^2/B^2)
  integer, protected, allocatable      :: grad_kappa_b(:)
  !> Variable index for (B . grad)B (curvature B drift)
  integer, protected, allocatable      :: b_dot_grad_b(:)
 
  ! ExB related drifts (vE = ExB/B^2)
  !> Variable index for curvature drift
  integer, protected, allocatable      :: ve_dot_grad_b(:)
  !> Variable index for polarization drift
  integer, protected, allocatable      :: b_dot_grad_ve(:)
  !> Variable index for polarization drift
  integer, protected, allocatable      :: ve_dot_grad_ve(:)
 
  public :: gca_init
  public :: gca_create_particles
  integer, parameter :: rk4=1, ark4=2
 
  ! Variables
  public :: bp, ep, grad_kappa_b, b_dot_grad_b
  public :: ve_dot_grad_b, b_dot_grad_ve, ve_dot_grad_ve
  public :: vp, jp, rhop
 
contains
 
  subroutine gca_init()
    use mod_global_parameters
    integer :: idir, nwx
 
    if (physics_type/='mhd') call mpistop("GCA particles need magnetic field!")
    if (ndir/=3) call mpistop("GCA particles need ndir=3!")
 
    nwx = 0
    allocate(bp(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      bp(idir) = nwx
    end do
    allocate(ep(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      ep(idir) = nwx
    end do
    allocate(grad_kappa_b(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      grad_kappa_b(idir) = nwx
    end do
    allocate(b_dot_grad_b(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      b_dot_grad_b(idir) = nwx
    end do
    allocate(ve_dot_grad_b(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      ve_dot_grad_b(idir) = nwx
    end do
    allocate(b_dot_grad_ve(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      b_dot_grad_ve(idir) = nwx
    end do
    allocate(ve_dot_grad_ve(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      ve_dot_grad_ve(idir) = nwx
    end do
    allocate(vp(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      vp(idir) = nwx
    end do
    allocate(jp(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      jp(idir) = nwx
    end do
    nwx = nwx + 1 ! density
    rhop = nwx
 
    ngridvars=nwx
 
    particles_fill_gridvars => gca_fill_gridvars
 
    if (associated(particles_define_additional_gridvars)) then
      call particles_define_additional_gridvars(ngridvars)
    end if
 
    select case(integrator_type_particles)
    case('RK4','Rk4','rk4')
      integrator = rk4
    case('ARK4','ARk4','Ark4','ark4')
      integrator = ark4
    case default
      integrator = ark4
    end select
 
    particles_integrate => gca_integrate_particles
 
  end subroutine gca_init
 
  subroutine gca_create_particles()
    ! initialise the particles
    use mod_global_parameters
    use mod_usr_methods, only: usr_create_particles, usr_update_payload, usr_check_particle
 
    double precision :: b(ndir), u(ndir), magmom
    double precision :: bnorm, lfac, vnorm, vperp, vpar
    double precision :: x(3, num_particles)
    double precision :: v(3, num_particles)
    double precision :: q(num_particles)
    double precision :: m(num_particles)
    double precision :: rrd(num_particles,ndir)
    double precision :: defpayload(ndefpayload)
    double precision :: usrpayload(nusrpayload)
    integer          :: igrid, ipe_particle
    integer          :: n, idir, nparticles_local
    logical          :: follow(num_particles), check
 
    if (mype==0) then
      if (.not. associated(usr_create_particles)) then
        ! Randomly distributed
        do idir=1,ndir
          do n = 1, num_particles
            rrd(n,idir) = rng%unif_01()
          end do
        end do
        do n=1, num_particles
          {^d&x(^d,n) = xprobmin^d + rrd(n+1,^d) * (xprobmax^d - xprobmin^d)\}
        end do
      else
        call usr_create_particles(num_particles, x, v, q, m, follow)
      end if
    end if
 
    call mpi_bcast(x,3*num_particles,mpi_double_precision,0,icomm,ierrmpi)
    call mpi_bcast(v,3*num_particles,mpi_double_precision,0,icomm,ierrmpi)
    call mpi_bcast(q,num_particles,mpi_double_precision,0,icomm,ierrmpi)
    call mpi_bcast(m,num_particles,mpi_double_precision,0,icomm,ierrmpi)
    call mpi_bcast(follow,num_particles,mpi_logical,0,icomm,ierrmpi)
 
    nparticles_local = 0
 
    ! first find ipe and igrid responsible for particle
    do n = 1, num_particles
      call find_particle_ipe(x(:, n),igrid,ipe_particle)
 
      particle(n)%igrid = igrid
      particle(n)%ipe   = ipe_particle
 
      if(ipe_particle == mype) then
        check = .true.
 
        ! Check for user-defined modifications or rejection conditions
        if (associated(usr_check_particle)) call usr_check_particle(igrid, x(:,n), v(:,n), q(n), m(n), follow(n), check)
        if (check) then
          call push_particle_into_particles_on_mype(n)
        else
          cycle
        end if
 
        nparticles_local = nparticles_local + 1
 
        call get_lfac_from_velocity(v(:, n), lfac)
 
        allocate(particle(n)%self)
        particle(n)%self%x      = x(:, n)
        particle(n)%self%q      = q(n)
        particle(n)%self%m      = m(n)
        particle(n)%self%follow = follow(n)
        particle(n)%self%index  = n
        particle(n)%self%time   = global_time
        particle(n)%self%dt     = 0.0d0
 
        call get_vec(bp, igrid, x(:, n), particle(n)%self%time, b)
 
        bnorm = norm2(b(:))
        vnorm = norm2(v(:, n))
        vpar  = sum(v(:, n) * b/bnorm)
        vperp = sqrt(vnorm**2 - vpar**2)
 
        ! The momentum vector u(1:3) is filled with the following components
 
        ! parallel momentum component (gamma v||)
        particle(n)%self%u(1) = lfac * vpar
 
        ! Mr: the conserved magnetic moment
        magmom = m(n) * (vperp * lfac)**2 / (2.0d0 * bnorm)
        particle(n)%self%u(2) = magmom
 
        ! Lorentz factor
        particle(n)%self%u(3) = lfac
 
        ! initialise payloads for GCA module
        allocate(particle(n)%payload(npayload))
        call gca_update_payload(igrid,x(:,n),particle(n)%self%u(:),q(n),m(n),defpayload,ndefpayload,0.d0)
        particle(n)%payload(1:ndefpayload) = defpayload
        if (associated(usr_update_payload)) then
          call usr_update_payload(igrid,x(:,n),particle(n)%self%u(:),q(n),m(n),usrpayload,nusrpayload,0.d0)
          particle(n)%payload(ndefpayload+1:npayload) = usrpayload
        end if
      end if
    end do
 
    call mpi_allreduce(nparticles_local,nparticles,1,mpi_integer,mpi_sum,icomm,ierrmpi)
 
  end subroutine gca_create_particles
 
  subroutine gca_fill_gridvars
    use mod_global_parameters
    use mod_usr_methods, only: usr_particle_fields
    use mod_geometry
 
    double precision, dimension(ixG^T,1:ndir) :: beta
    double precision, dimension(ixG^T,1:nw)   :: w,wold
    double precision                          :: current(ixG^T,7-2*ndir:3)
    double precision, dimension(ixG^T,1:ndir) :: vE, bhat
    double precision, dimension(ixG^T)        :: kappa, kappa_B, absB, tmp
    integer                                   :: igrid, iigrid, idir, idim
    integer                                   :: idirmin
 
    call fill_gridvars_default()
 
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      w(ixg^t,1:nw) = ps(igrid)%w(ixg^t,1:nw)
      call phys_to_primitive(ixg^ll,ixg^ll,w,ps(igrid)%x)
 
      ! grad(kappa B)
      absb(ixg^t) = sqrt(sum(gridvars(igrid)%w(ixg^t,bp(:))**2,dim=ndim+1))
      ve(ixg^t,1) = gridvars(igrid)%w(ixg^t,ep(2)) * gridvars(igrid)%w(ixg^t,bp(3)) &
            - gridvars(igrid)%w(ixg^t,ep(3)) * gridvars(igrid)%w(ixg^t,bp(2))
      ve(ixg^t,2) = gridvars(igrid)%w(ixg^t,ep(3)) * gridvars(igrid)%w(ixg^t,bp(1)) &
            - gridvars(igrid)%w(ixg^t,ep(1)) * gridvars(igrid)%w(ixg^t,bp(3))
      ve(ixg^t,3) = gridvars(igrid)%w(ixg^t,ep(1)) * gridvars(igrid)%w(ixg^t,bp(2)) &
            - gridvars(igrid)%w(ixg^t,ep(2)) * gridvars(igrid)%w(ixg^t,bp(1))
      do idir=1,ndir
        where (absb(ixg^t) .gt. 0.d0)
          ve(ixg^t,idir) = ve(ixg^t,idir) / absb(ixg^t)**2
        elsewhere
          ve(ixg^t,idir) = 0.d0
        end where
      end do
      if (any(sum(ve(ixg^t,:)**2,dim=ndim+1) .ge. c_norm**2) .and. relativistic) then
        call mpistop("GCA FILL GRIDVARS: vE>c! ABORTING...")
      end if
      if (any(ve .ne. ve)) then
        call mpistop("GCA FILL GRIDVARS: NaNs IN vE! ABORTING...")
      end if
 
      if (relativistic) then
        kappa(ixg^t) = 1.d0/sqrt(1.0d0 - sum(ve(ixg^t,:)**2,dim=ndim+1)/c_norm**2)
      else
        kappa(ixg^t) = 1.d0
      end if
      kappa_b(ixg^t) = absb(ixg^t) / kappa(ixg^t)
 
      if (any(kappa_b .ne. kappa_b)) then
        call mpistop("GCA FILL GRIDVARS: NaNs IN kappa_B! ABORTING...")
      end if
 
      tmp=0.d0
      do idim=1,ndim
        call gradient(kappa_b,ixg^ll,ixg^ll^lsub1,idim,tmp)
        gridvars(igrid)%w(ixg^t,grad_kappa_b(idim)) = tmp(ixg^t)
      end do
 
      ! bhat
      do idir=1,ndir
        where (absb(ixg^t) .gt. 0.d0)
          bhat(ixg^t,idir) = gridvars(igrid)%w(ixg^t,bp(idir)) / absb(ixg^t)
        elsewhere
          bhat(ixg^t,idir) = 0.d0
        end where
      end do
      if (any(bhat .ne. bhat)) then
        call mpistop("GCA FILL GRIDVARS: NaNs IN bhat! ABORTING...")
      end if
 
      do idir=1,ndir
        ! (b dot grad) b and the other directional derivatives
        do idim=1,ndim
          call gradient(bhat(ixg^t,idir),ixg^ll,ixg^ll^lsub1,idim,tmp)
          gridvars(igrid)%w(ixg^t,b_dot_grad_b(idir)) = gridvars(igrid)%w(ixg^t,b_dot_grad_b(idir)) &
               + bhat(ixg^t,idim) * tmp(ixg^t)
          gridvars(igrid)%w(ixg^t,ve_dot_grad_b(idir)) = gridvars(igrid)%w(ixg^t,ve_dot_grad_b(idir)) &
               + ve(ixg^t,idim) * tmp(ixg^t)
          call gradient(ve(ixg^t,idir),ixg^ll,ixg^ll^lsub1,idim,tmp)
          gridvars(igrid)%w(ixg^t,b_dot_grad_ve(idir)) = gridvars(igrid)%w(ixg^t,b_dot_grad_ve(idir)) &
               + bhat(ixg^t,idim) * tmp(ixg^t)
          gridvars(igrid)%w(ixg^t,ve_dot_grad_ve(idir)) = gridvars(igrid)%w(ixg^t,ve_dot_grad_ve(idir)) &
               + ve(ixg^t,idim) * tmp(ixg^t)
        end do
      end do
 
    end do
 
  end subroutine gca_fill_gridvars
 
  subroutine gca_integrate_particles(end_time)
    use mod_odeint
    use mod_global_parameters
    use mod_usr_methods, only: usr_create_particles, usr_update_payload
    double precision, intent(in)        :: end_time
 
    double precision                    :: lfac, absS
    double precision                    :: defpayload(ndefpayload)
    double precision                    :: usrpayload(nusrpayload)
    double precision                    :: dt_p, tloc, y(ndir+2),dydt(ndir+2),ytmp(ndir+2), euler_cfl, int_factor
    double precision                    :: tk, k1(ndir+2),k2(ndir+2),k3(ndir+2),k4(ndir+2)
    double precision, dimension(1:ndir) :: x, vE, e, b, bhat, x_new, vfluid, current
    double precision, dimension(1:ndir) :: drift1, drift2
    double precision, dimension(1:ndir) :: drift3, drift4, drift5, drift6, drift7
    double precision, dimension(1:ndir) :: bdotgradb, vEdotgradb, gradkappaB
    double precision, dimension(1:ndir) :: bdotgradvE, vEdotgradvE
    double precision, dimension(1:ndir) :: gradBdrift, reldrift, bdotgradbdrift
    double precision, dimension(1:ndir) :: vEdotgradbdrift, bdotgradvEdrift
    double precision, dimension(1:ndir) :: vEdotgradvEdrift
    double precision                    :: kappa, Mr, upar, m, absb, gamma, q, mompar, vpar, vEabs
    double precision                    :: gradBdrift_abs, reldrift_abs, epar
    double precision                    :: bdotgradbdrift_abs, vEdotgradbdrift_abs
    double precision                    :: bdotgradvEdrift_abs, vEdotgradvEdrift_abs
    double precision                    :: momentumpar1, momentumpar2, momentumpar3, momentumpar4
    ! Precision of time-integration:
    double precision,parameter          :: eps=1.0d-6
    ! for odeint:
    double precision                    :: h1, hmin, h_old
    integer                             :: nok, nbad, ic1^D, ic2^D, ierror, nvar
    integer                             :: ipart, iipart, seed, ic^D,igrid_particle, ipe_particle, ipe_working
    logical                             :: int_choice
    logical                             :: BC_applied
 
    nvar=ndir+2
 
    do iipart=1,nparticles_active_on_mype
      ipart                   = particles_active_on_mype(iipart)
      int_choice              = .false.
      dt_p                    = gca_get_particle_dt(particle(ipart), end_time)
      particle(ipart)%self%dt = dt_p
 
      igrid_working           = particle(ipart)%igrid
      ipart_working           = particle(ipart)%self%index
      tloc                    = particle(ipart)%self%time
      x(1:ndir)               = particle(ipart)%self%x(1:ndir)
 
      select case (integrator)
      case (rk4)
        ! Runge-Kutta order 4
        tk = tloc
        y(1:ndir) = x(1:ndir) ! position of guiding center
        y(ndir+1) = particle(ipart)%self%u(1) ! parallel momentum component (gamma v||)
        y(ndir+2) = particle(ipart)%self%u(2) ! conserved magnetic moment Mr
        call derivs_gca_rk(tk,y,k1)
        tk = tloc + dt_p/2.d0
        y(1:ndir) = x + dt_p*k1(1:ndir)/2.d0
        y(ndir+1) = particle(ipart)%self%u(1) + dt_p*k1(ndir+1)/2.d0
        call derivs_gca_rk(tk,y,k2)
        tk = tloc + dt_p/2.d0
        y(1:ndir) = x + dt_p*k2(1:ndir)/2.d0
        y(ndir+1) = particle(ipart)%self%u(1) + dt_p*k2(ndir+1)/2.d0
        call derivs_gca_rk(tk,y,k3)
        tk = tloc + dt_p
        y(1:ndir) = x + dt_p*k3(1:ndir)
        y(ndir+1) = particle(ipart)%self%u(1) + dt_p*k3(ndir+1)
        call derivs_gca_rk(tk,y,k4)
        y(1:ndir) = x(1:ndir) ! position of guiding center
        y(ndir+1) = particle(ipart)%self%u(1) ! parallel momentum component (gamma v||)
        y = y + dt_p/6.d0*(k1 + 2.d0*k2 + 2.d0*k3 + k4)
        particle(ipart)%self%x(1:ndir) = y(1:ndir)
        particle(ipart)%self%u(1)      = y(ndir+1)
 
      case (ark4)
        ! Adaptive stepwidth RK4:
        ! initial solution vector:
        y(1:ndir) = x(1:ndir) ! position of guiding center
        y(ndir+1) = particle(ipart)%self%u(1) ! parallel momentum component (gamma v||)
        y(ndir+2) = particle(ipart)%self%u(2) ! conserved magnetic moment Mr
      ! y(ndir+3) = particle(ipart)%self%u(3) ! Lorentz factor of particle
  
        ! we temporarily save the solution vector, to replace the one from the euler
        ! timestep after euler integration
        ytmp=y
  
        call derivs_gca(particle(ipart)%self%time,y,dydt)
  
        ! make an Euler step with the proposed timestep:
        ! factor to ensure we capture all particles near the internal ghost cells.
        ! Can be adjusted during a run, after an interpolation error.
        euler_cfl=2.5d0
  
        ! new solution vector:
        y(1:ndir+2) = y(1:ndir+2) + euler_cfl * dt_p * dydt(1:ndir+2)
        particle(ipart)%self%x(1:ndir) = y(1:ndir) ! position of guiding center
        particle(ipart)%self%u(1)      = y(ndir+1) ! parallel momentum component(gamma v||)
        particle(ipart)%self%u(2)      = y(ndir+2) ! conserved magnetic moment
  
        ! check if the particle is in the internal ghost cells
        int_factor =1.0d0
  
        if(.not. particle_in_igrid(ipart_working,igrid_working)) then
          ! if particle is not in the grid an euler timestep is taken instead of a RK4
          ! timestep. Then based on that we do an interpolation and check how much further
          ! the timestep for the RK4 has to be restricted.
          ! factor to make integration more accurate for particles near the internal
          ! ghost cells. This factor can be changed during integration after an
          ! interpolation error. But one should be careful with timesteps for i/o
  
          ! flat interpolation:
          {ic^d = int((y(^d)-rnode(rpxmin^d_,igrid_working))/rnode(rpdx^d_,igrid_working)) + 1 + nghostcells\}
  
          ! linear interpolation:
          {
          if (ps(igrid_working)%x({ic^dd},^d) .lt. y(^d)) then
             ic1^d = ic^d
          else
             ic1^d = ic^d -1
          end if
          ic2^d = ic1^d + 1
          \}
  
          int_factor =0.5d0
  
          {^d&
          if (ic1^d .le. ixglo^d-2 .or. ic2^d .ge. ixghi^d+2) then
            int_factor = 0.05d0
          end if
          \}
  
          {^d&
          if (ic1^d .eq. ixglo^d-1 .or. ic2^d .eq. ixghi^d+1) then
            int_factor = 0.1d0
          end if
          \}
  
          dt_p=int_factor*dt_p
        end if
  
        ! replace the solution vector with the original as it was before the Euler timestep
        y(1:ndir+2) = ytmp(1:ndir+2)
  
        particle(ipart)%self%x(1:ndir) = ytmp(1:ndir) ! position of guiding center
        particle(ipart)%self%u(1)      = ytmp(ndir+1) ! parallel momentum component (gamma v||)
        particle(ipart)%self%u(2)      = ytmp(ndir+2) ! conserved magnetic moment
  
        ! specify a minimum step hmin. If the timestep reaches this minimum, multiply by
        ! a factor 100 to make sure the RK integration doesn't crash
        h1 = dt_p/2.0d0; hmin=1.0d-9; h_old=dt_p/2.0d0
  
        if(h1 .lt. hmin)then
          h1=hmin
          dt_p=2.0d0*h1
        endif
  
        if (any(y .ne. y)) then
          call mpistop("NaNs DETECTED IN GCA_INTEGRATE BEFORE ODEINT CALL! ABORTING...")
        end if
  
        ! RK4 integration with adaptive stepwidth
        call odeint(y,nvar,tloc,tloc+dt_p,eps,h1,hmin,nok,nbad,derivs_gca_rk,rkqs,ierror)
  
        if (ierror /= 0) then
           print *, "odeint returned error code", ierror
           print *, "1 means hmin too small, 2 means MAXSTP exceeded"
           print *, "Having a problem with particle", iipart
        end if
  
        if (any(y .ne. y)) then
          call mpistop("NaNs DETECTED IN GCA_INTEGRATE AFTER ODEINT CALL! ABORTING...")
        end if
  
        ! original RK integration without interpolation in ghost cells
  !       call odeint(y,nvar,tloc,tloc+dt_p,eps,h1,hmin,nok,nbad,derivs_gca,rkqs)
  
        ! final solution vector after rk integration
        particle(ipart)%self%x(1:ndir) = y(1:ndir)
        particle(ipart)%self%u(1)      = y(ndir+1)
        particle(ipart)%self%u(2)      = y(ndir+2)
        !particle(ipart)%self%u(3)      = y(ndir+3)
      end select
 
      ! now calculate other quantities, mean Lorentz factor, drifts, perpendicular velocity:
      call get_bfield(igrid_working, y(1:ndir), tloc+dt_p, b)
      call get_efield(igrid_working, y(1:ndir), tloc+dt_p, e)
!      call get_vec(bp, igrid_working,y(1:ndir),tloc+dt_p,b)
!      call get_vec(vp, igrid_working,y(1:ndir),tloc+dt_p,vfluid)
!      call get_vec(jp, igrid_working,y(1:ndir),tloc+dt_p,current)
!      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)
 
      absb         = sqrt(sum(b(:)**2))
      bhat(1:ndir) = b(1:ndir) / absb
 
      epar         = sum(e(:)*bhat(:))
 
      call cross(e,bhat,ve)
 
      ve(1:ndir) = ve(1:ndir) / absb
      veabs = sqrt(sum(ve(:)**2))
      if (relativistic) then
        kappa = 1.d0/sqrt(1.0d0 - sum(ve(:)**2)/c_norm**2)
      else
        kappa = 1.d0
      end if
      mr = y(ndir+2); upar = y(ndir+1); m=particle(ipart)%self%m; q=particle(ipart)%self%q
      if (relativistic) then
        gamma = sqrt(1.0d0+upar**2/c_norm**2+2.0d0*mr*absb/m/c_norm**2)*kappa
      else
        gamma = 1.d0
      end if
 
      particle(ipart)%self%u(3) = gamma
 
      ! Time update
      particle(ipart)%self%time = particle(ipart)%self%time + dt_p
 
      ! Update payload
      call gca_update_payload(igrid_working,&
             particle(ipart)%self%x,particle(ipart)%self%u,q,m,defpayload,ndefpayload,particle(ipart)%self%time)
      particle(ipart)%payload(1:ndefpayload) = defpayload
      if (associated(usr_update_payload)) then
        call usr_update_payload(igrid_working,&
             particle(ipart)%self%x,particle(ipart)%self%u,q,m,usrpayload,nusrpayload,particle(ipart)%self%time)
        particle(ipart)%payload(ndefpayload+1:npayload) = usrpayload
      end if
 
    end do
 
  end subroutine gca_integrate_particles
 
  subroutine derivs_gca_rk(t_s,y,dydt)
    use mod_global_parameters
 
    double precision                :: t_s, y(ndir+2)
    double precision                :: dydt(ndir+2)
 
    double precision,dimension(ndir):: vE, b, e, x, bhat, bdotgradb, vEdotgradb, gradkappaB, vfluid, current
    double precision,dimension(ndir):: bdotgradvE, vEdotgradvE, u, utmp1, utmp2, utmp3
    double precision                :: upar, Mr, gamma, absb, q, m, epar, kappa
    integer                         :: ic^D
 
    ! Here the terms in the guiding centre equations of motion are interpolated for
    ! the RK integration. The interpolation is also done in the ghost cells such
    ! that the RK integration does not give an error
 
    q = particle(ipart_working)%self%q
    m = particle(ipart_working)%self%m
 
    x(1:ndir) = y(1:ndir)
    upar      = y(ndir+1) ! gamma v||
    mr        = y(ndir+2)
    !gamma     = y(ndir+3)
 
    if (any(x .ne. x)) then
      write(*,*) "ERROR IN DERIVS_GCA_RK: NaNs IN X OR Y!"
      write(*,*) "x",x
      write(*,*) "y",y(ndir+1:ndir+2)
      call mpistop("ABORTING...")
    end if
 
    call get_bfield(igrid_working, x, t_s, b)
    call get_efield(igrid_working, x, t_s, e)
!    call get_vec(bp, igrid_working,x,t_s,b)
!    call get_vec(vp, igrid_working,x,t_s,vfluid)
!    call get_vec(jp, igrid_working,x,t_s,current)
!    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)
    call get_vec(b_dot_grad_b, igrid_working,x,t_s,bdotgradb)
    call get_vec(ve_dot_grad_b, igrid_working,x,t_s,vedotgradb)
    call get_vec(grad_kappa_b, igrid_working,x,t_s,gradkappab)
    call get_vec(b_dot_grad_ve, igrid_working,x,t_s,bdotgradve)
    call get_vec(ve_dot_grad_ve, igrid_working,x,t_s,vedotgradve)
 
    if (any(b .ne. b) .or. any(e .ne. e) &
        .or. any(bdotgradb .ne. bdotgradb) .or. any(vedotgradb .ne. vedotgradb) &
        .or. any(gradkappab .ne. gradkappab) .or. any(bdotgradve .ne. bdotgradve) &
        .or. any(vedotgradve .ne. vedotgradve)) then
      write(*,*) "ERROR IN DERIVS_GCA_RK: NaNs IN FIELD QUANTITIES!"
      write(*,*) "b",b
      write(*,*) "e",e
      write(*,*) "bdotgradb",bdotgradb
      write(*,*) "vEdotgradb",vedotgradb
      write(*,*) "gradkappaB",gradkappab
      write(*,*) "bdotgradvE",bdotgradve
      write(*,*) "vEdotgradvE",vedotgradve
      call mpistop("ABORTING...")
    end if
 
    absb  = sqrt(sum(b(:)**2))
    if (absb .gt. 0.d0) then
      bhat(1:ndir) = b(1:ndir) / absb
    else
      bhat = 0.d0
    end if
    epar = sum(e(:)*bhat(:))
 
    call cross(e,bhat,ve)
    if (absb .gt. 0.d0) then
      ve(1:ndir) = ve(1:ndir) / absb
    else
      ve(1:ndir) = 0.d0
    end if
 
    if (relativistic) then
      kappa = 1.d0/sqrt(1.0d0 - sum(ve(:)**2)/c_norm**2)
      gamma = sqrt(1.0d0+upar**2/c_norm**2+2.0d0*mr*absb/m/c_norm**2)*kappa
    else
      kappa = 1.d0
      gamma = 1.d0
    end if
 
    if (absb .gt. 0.d0) then
      utmp1(1:ndir) = bhat(1:ndir)/(absb/kappa**2)
    else
      utmp1 = 0.d0
    end if
    utmp2(1:ndir) = mr/(gamma*q)*gradkappab(1:ndir) &
         + m/q* (upar**2/gamma*bdotgradb(1:ndir) + upar*vedotgradb(1:ndir) &
                 + upar*bdotgradve(1:ndir) + gamma*vedotgradve(1:ndir))
    if (relativistic) then
      utmp2(1:ndir) = utmp2(1:ndir) + upar*epar/(gamma)*ve(1:ndir)
    end if
 
    call cross(utmp1,utmp2,utmp3)
    u(1:ndir) = ve(1:ndir) + utmp3(1:ndir)
 
    ! done assembling the terms, now write rhs:
    dydt(1:ndir) = ( u(1:ndir) + upar/gamma * bhat(1:ndir) )
    dydt(ndir+1) = q/m*epar - mr/(m*gamma) * sum(bhat(:)*gradkappab(:)) &
                   + sum(ve(:)*(upar*bdotgradb(:)+gamma*vedotgradb(:)))
    dydt(ndir+2) = 0.0d0 ! magnetic moment is conserved
 
  end subroutine derivs_gca_rk
 
  subroutine derivs_gca(t_s,y,dydt)
    use mod_global_parameters
 
    double precision                :: t_s, y(ndir+2)
    double precision                :: dydt(ndir+2)
 
    double precision,dimension(ndir):: vE, b, e, x, bhat, bdotgradb, vEdotgradb, gradkappaB, vfluid, current
    double precision,dimension(ndir):: bdotgradvE, vEdotgradvE, u, utmp1, utmp2, utmp3
    double precision                :: upar, Mr, gamma, absb, q, m, epar, kappa
 
    ! Here the normal interpolation is done for the terms in the GCA equations of motion
 
    q = particle(ipart_working)%self%q
    m = particle(ipart_working)%self%m
 
    x(1:ndir) = y(1:ndir)
    upar      = y(ndir+1) ! gamma v||
    mr        = y(ndir+2)
    !gamma     = y(ndir+3)
 
    if (any(x .ne. x)) then
      call mpistop("ERROR IN DERIVS_GCA: NaNs IN X! ABORTING...")
    end if
 
    call get_bfield(igrid_working, x, t_s, b)
    call get_efield(igrid_working, x, t_s, e)
!    call get_vec(bp, igrid_working,x,t_s,b)
!    call get_vec(vp, igrid_working,x,t_s,vfluid)
!    call get_vec(jp, igrid_working,x,t_s,current)
!    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)
    call get_vec(b_dot_grad_b, igrid_working,x,t_s,bdotgradb)
    call get_vec(ve_dot_grad_b, igrid_working,x,t_s,vedotgradb)
    call get_vec(grad_kappa_b, igrid_working,x,t_s,gradkappab)
    call get_vec(b_dot_grad_ve, igrid_working,x,t_s,bdotgradve)
    call get_vec(ve_dot_grad_ve, igrid_working,x,t_s,vedotgradve)
 
    absb         = sqrt(sum(b(:)**2))
    if (absb .gt. 0.d0) then
      bhat(1:ndir) = b(1:ndir) / absb
    else
      bhat = 0.d0
    end if
 
    epar         = sum(e(:)*bhat(:))
    call cross(e,bhat,ve)
    if (absb .gt. 0.d0) ve(1:ndir)   = ve(1:ndir) / absb
 
    if (relativistic) then
      kappa = 1.d0/sqrt(1.0d0 - sum(ve(:)**2)/c_norm**2)
      gamma = sqrt(1.0d0+upar**2/c_norm**2+2.0d0*mr*absb/m/c_norm**2)*kappa
    else
      kappa = 1.d0
      gamma = 1.d0
    end if
    if (absb .gt. 0.d0) then
      utmp1(1:ndir) = bhat(1:ndir)/(absb/kappa**2)
    else
      utmp1 = 0.d0
    end if
    utmp2(1:ndir) = mr/(gamma*q)*gradkappab(1:ndir) &
         + m/q* (upar**2/gamma*bdotgradb(1:ndir) + upar*vedotgradb(1:ndir) &
                 + upar*bdotgradve(1:ndir) + gamma*vedotgradve(1:ndir))
    if (relativistic) then
      utmp2(1:ndir) = utmp2(1:ndir) + upar*epar/(gamma)*ve(1:ndir)
    end if
 
    call cross(utmp1,utmp2,utmp3)
    u(1:ndir) = ve(1:ndir) + utmp3(1:ndir)
 
    ! done assembling the terms, now write rhs:
    dydt(1:ndir) = ( u(1:ndir) + upar/gamma * bhat(1:ndir) )
    dydt(ndir+1) = q/m*epar - mr/(m*gamma) * sum(bhat(:)*gradkappab(:)) &
                   + sum(ve(:)*(upar*bdotgradb(:)+gamma*vedotgradb(:)))
    dydt(ndir+2) = 0.0d0 ! magnetic moment is conserved
 
  end subroutine derivs_gca
 
  !> Update payload subroutine
  subroutine gca_update_payload(igrid,xpart,upart,qpart,mpart,mypayload,mynpayload,particle_time)
    use mod_global_parameters
    integer, intent(in)           :: igrid,mynpayload
    double precision, intent(in)  :: xpart(1:ndir),upart(1:ndir),qpart,mpart,particle_time
    double precision, intent(out) :: mypayload(mynpayload)
    double precision, dimension(1:ndir) :: vE, e, b, bhat, vfluid, current
    double precision, dimension(1:ndir) :: drift1, drift2
    double precision, dimension(1:ndir) :: drift3, drift4, drift5, drift6, drift7
    double precision, dimension(1:ndir) :: bdotgradb, vEdotgradb, gradkappaB
    double precision, dimension(1:ndir) :: bdotgradvE, vEdotgradvE
    double precision, dimension(1:ndir) :: gradBdrift, reldrift, bdotgradbdrift
    double precision, dimension(1:ndir) :: vEdotgradbdrift, bdotgradvEdrift
    double precision, dimension(1:ndir) :: vEdotgradvEdrift
    double precision                    :: kappa, upar, absb, gamma, vpar, vEabs
    double precision                    :: gradBdrift_abs, reldrift_abs, epar
    double precision                    :: bdotgradbdrift_abs, vEdotgradbdrift_abs
    double precision                    :: bdotgradvEdrift_abs, vEdotgradvEdrift_abs
    double precision                    :: momentumpar1, momentumpar2, momentumpar3, momentumpar4
 
    call get_bfield(igrid, xpart(1:ndir), particle_time, b)
    call get_efield(igrid, xpart(1:ndir), particle_time, e)
!    call get_vec(bp, igrid,xpart(1:ndir),particle_time,b)
!    call get_vec(vp, igrid,xpart(1:ndir),particle_time,vfluid)
!    call get_vec(jp, igrid,xpart(1:ndir),particle_time,current)
!    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)
 
    absb         = sqrt(sum(b(:)**2))
    bhat(1:ndir) = b(1:ndir) / absb
    epar         = sum(e(:)*bhat(:))
    call cross(e,bhat,ve)
    ve(1:ndir)   = ve(1:ndir) / absb
    veabs = sqrt(sum(ve(:)**2))
    if (relativistic) then
      kappa = 1.d0/sqrt(1.0d0 - sum(ve(:)**2)/c_norm**2)
    else
      kappa = 1.d0
    end if
    vpar = upart(1)/upart(3)
    upar = upart(1)
 
    call get_vec(b_dot_grad_b, igrid,xpart(1:ndir),particle_time,bdotgradb)
    call get_vec(ve_dot_grad_b, igrid,xpart(1:ndir),particle_time,vedotgradb)
    call get_vec(grad_kappa_b, igrid,xpart(1:ndir),particle_time,gradkappab)
    call get_vec(b_dot_grad_ve, igrid,xpart(1:ndir),particle_time,bdotgradve)
    call get_vec(ve_dot_grad_ve, igrid,xpart(1:ndir),particle_time,vedotgradve)
 
    drift1(1:ndir) = bhat(1:ndir)/(absb/kappa**2)
    drift2(1:ndir) = upart(2)/(upart(3)*qpart)*gradkappab(1:ndir)
 
    call cross(drift1,drift2,gradbdrift)
    gradbdrift_abs = sqrt(sum(gradbdrift(:)**2))
 
    drift3(1:ndir) = upar*epar/upart(3)*ve(1:ndir)
    call cross(drift1,drift3,reldrift)
    reldrift_abs = sqrt(sum(reldrift(:)**2))
 
    drift4(1:ndir) = mpart/qpart* ( upar**2/upart(3)*bdotgradb(1:ndir))
    call cross(drift1,drift4,bdotgradbdrift)
    bdotgradbdrift_abs = sqrt(sum(bdotgradbdrift(:)**2))
 
    drift5(1:ndir) = mpart/qpart* ( upar*vedotgradb(1:ndir))
    call cross(drift1,drift5,vedotgradbdrift)
    vedotgradbdrift_abs = sqrt(sum(vedotgradbdrift(:)**2))
 
    drift6(1:ndir) = mpart/qpart* ( upar*bdotgradve(1:ndir))
    call cross(drift1,drift6,bdotgradvedrift)
    bdotgradvedrift_abs = sqrt(sum(bdotgradvedrift(:)**2))
 
    drift7(1:ndir) = mpart/qpart* (upart(3)*vedotgradve(1:ndir))
    call cross(drift1,drift7,vedotgradvedrift)
    vedotgradvedrift_abs = sqrt(sum(vedotgradvedrift(:)**2))
 
    momentumpar1 = qpart/mpart*epar
    momentumpar2 = -(upart(2)/mpart/upart(3))*sum(bhat(:)*gradkappab(:))
    momentumpar3 = upar*sum(ve(:)*bdotgradb(:))
    momentumpar4 = upart(3)*sum(ve(:)*vedotgradb(:))
 
    ! Payload update
    if (mynpayload > 0) then
      ! current gyroradius
      mypayload(1) = sqrt(2.0d0*mpart*upart(2)*absb)/abs(qpart*absb)
    end if
    if (mynpayload > 1) then
      ! pitch angle
      mypayload(2) = atan2(sqrt((2.0d0*upart(2)*absb)/(mpart*upart(3)**2)),vpar)
    end if
    if (mynpayload > 2) then
      ! particle v_perp
      mypayload(3) = sqrt((2.0d0*upart(2)*absb)/(mpart*upart(3)**2))
    end if
    if (mynpayload > 3) then
      ! particle parallel momentum term 1
      mypayload(4) = momentumpar1
    end if
    if (mynpayload > 4) then
      ! particle parallel momentum term 2
      mypayload(5) = momentumpar2
    end if
    if (mynpayload > 5) then
      ! particle parallel momentum term 3
      mypayload(6) = momentumpar3
    end if
    if (mynpayload > 6) then
      ! particle parallel momentum term 4
      mypayload(7) = momentumpar4
    end if
    if (mynpayload > 7) then
      ! particle ExB drift
      mypayload(8) = veabs
    end if
    if (mynpayload > 8) then
      ! relativistic drift
      mypayload(9) = gradbdrift_abs
    end if
    if (mynpayload > 9) then
      ! gradB drift
      mypayload(10) = reldrift_abs
    end if
    if (mynpayload > 10) then
      ! bdotgradb drift
      mypayload(11) = bdotgradbdrift_abs
    end if
    if (mynpayload > 11) then
      ! vEdotgradb drift
      mypayload(12) = vedotgradbdrift_abs
    end if
    if (mynpayload > 12) then
      ! bdotgradvE drift
      mypayload(13) = bdotgradvedrift_abs
    end if
    if (mynpayload > 13) then
      ! vEdotgradvE drift
      mypayload(14) = vedotgradvedrift_abs
    end if
 
  end subroutine gca_update_payload
 
  function gca_get_particle_dt(partp, end_time) result(dt_p)
    use mod_odeint
    use mod_global_parameters
    type(particle_ptr), intent(in) :: partp
    double precision, intent(in)   :: end_time
    double precision               :: dt_p
 
    double precision            :: tout, dt_particles_mype, dt_cfl0, dt_cfl1, dt_a
    double precision            :: dxmin, vp, a, gammap
    double precision            :: v(ndir), y(ndir+2),ytmp(ndir+2), dydt(ndir+2), v0(ndir), v1(ndir), dydt1(ndir+2)
    double precision            :: ap0, ap1, dt_cfl_ap0, dt_cfl_ap1, dt_cfl_ap
    double precision            :: dt_euler, dt_tmp
    ! make these particle cfl conditions more restrictive if you are interpolating out of the grid
    double precision            :: cfl, uparcfl
    double precision, parameter :: uparmin=1.0d-6*const_c
    integer                     :: ipart, iipart, nout, ic^d, igrid_particle, ipe_particle, ipe
    logical                     :: bc_applied
 
    if (const_dt_particles > 0) then
      dt_p = const_dt_particles
      return
    end if
 
    cfl = particles_cfl
    uparcfl = particles_cfl
 
    igrid_working = partp%igrid
    ipart_working = partp%self%index
    dt_tmp = (end_time - partp%self%time)
    if(dt_tmp .le. 0.0d0) dt_tmp = smalldouble
    ! make sure we step only one cell at a time, first check CFL at current location
    ! then we make an Euler step to the new location and check the new CFL
    ! we simply take the minimum of the two timesteps.
    ! added safety factor cfl:
    dxmin  = min({rnode(rpdx^d_,igrid_working)},bigdouble)*cfl
    ! initial solution vector:
    y(1:ndir) = partp%self%x(1:ndir) ! position of guiding center
    y(ndir+1) = partp%self%u(1) ! parallel momentum component (gamma v||)
    y(ndir+2) = partp%self%u(2) ! conserved magnetic moment
    ytmp=y
    !y(ndir+3) = partp%self%u(3) ! Lorentz factor of guiding centre
 
    call derivs_gca(partp%self%time,y,dydt)
    v0(1:ndir) = dydt(1:ndir)
    ap0        = dydt(ndir+1)
 
    ! guiding center velocity:
    v(1:ndir) = abs(dydt(1:ndir))
    vp = sqrt(sum(v(:)**2))
 
    dt_cfl0    = dxmin / max(vp, smalldouble)
    dt_cfl_ap0 = uparcfl * abs(max(abs(y(ndir+1)),uparmin) / max(abs(ap0), smalldouble))
    !dt_cfl_ap0 = min(dt_cfl_ap0, uparcfl * sqrt(abs(unit_length*dxmin/(ap0+smalldouble))) )
 
    ! make an Euler step with the proposed timestep:
    ! new solution vector:
    dt_euler = min(dt_tmp,dt_cfl0,dt_cfl_ap0)
    y(1:ndir+2) = y(1:ndir+2) + dt_euler * dydt(1:ndir+2)
 
    partp%self%x(1:ndir) = y(1:ndir) ! position of guiding center
    partp%self%u(1)      = y(ndir+1) ! parallel momentum component (gamma v||)
    partp%self%u(2)      = y(ndir+2) ! conserved magnetic moment
 
    ! first check if the particle is outside the physical domain or in the ghost cells
    if(.not. particle_in_igrid(ipart_working,igrid_working)) then
      y(1:ndir+2) = ytmp(1:ndir+2)
    end if
 
    call derivs_gca_rk(partp%self%time+dt_euler,y,dydt)
    !call derivs_gca(partp%self%time+dt_euler,y,dydt)
 
    v1(1:ndir) = dydt(1:ndir)
    ap1        = dydt(ndir+1)
 
    ! guiding center velocity:
    v(1:ndir) = abs(dydt(1:ndir))
    vp = sqrt(sum(v(:)**2))
 
    dt_cfl1    = dxmin / max(vp, smalldouble)
    dt_cfl_ap1 = uparcfl * abs(max(abs(y(ndir+1)),uparmin) / max(abs(ap1), smalldouble))
    !dt_cfl_ap1 = min(dt_cfl_ap1, uparcfl * sqrt(abs(unit_length*dxmin/(ap1+smalldouble))) )
 
    dt_tmp = min(dt_euler, dt_cfl1, dt_cfl_ap1)
 
    partp%self%x(1:ndir) = ytmp(1:ndir) ! position of guiding center
    partp%self%u(1)      = ytmp(ndir+1) ! parallel momentum component (gamma v||)
    partp%self%u(2)      = ytmp(ndir+2) ! conserved magnetic moment
    !dt_tmp = min(dt_cfl1, dt_cfl_ap1)
 
    ! time step due to parallel acceleration:
    ! The standard thing, dt=sqrt(dx/a) where we compute a from d(gamma v||)/dt and d(gamma)/dt
    ! dt_ap = sqrt(abs(dxmin*unit_length*y(ndir+3)/( dydt(ndir+1) - y(ndir+1)/y(ndir+3)*dydt(ndir+3) ) ) )
    ! vp = sqrt(sum(v(1:ndir)**))
    ! gammap = sqrt(1.0d0/(1.0d0-(vp/const_c)**2))
    ! ap = const_c**2/vp*gammap**(-3)*dydt(ndir+3)
    ! dt_ap = sqrt(dxmin*unit_length/ap)
 
    !dt_a = bigdouble
    !if (dt_euler .gt. smalldouble) then
    !   a = sqrt(sum((v1(1:ndir)-v0(1:ndir))**2))/dt_euler
    !   dt_a = min(sqrt(dxmin/a),bigdouble)
    !end if
 
    !dt_p = min(dt_tmp , dt_a)
    dt_p = dt_tmp
 
    ! Make sure we do not advance beyond end_time
    call limit_dt_endtime(end_time - partp%self%time, dt_p)
 
  end function gca_get_particle_dt
 
end module mod_particle_gca