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, ic11,ic12,ic13, ic21,&
ic22,ic23, ierror, nvar
integer :: ipart, iipart, seed, ic1,ic2,ic3,&
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:
ic1 = int((y(1)-rnode(rpxmin1_,igrid_working))/rnode(rpdx1_,&
igrid_working)) + 1 + nghostcells
ic2 = int((y(2)-rnode(rpxmin2_,igrid_working))/rnode(rpdx2_,&
igrid_working)) + 1 + nghostcells
ic3 = int((y(3)-rnode(rpxmin3_,igrid_working))/rnode(rpdx3_,&
igrid_working)) + 1 + nghostcells
! linear interpolation:
if (ps(igrid_working)%x(ic1,ic2,ic3,1) .lt. y(1)) then
ic11 = ic1
else
ic11 = ic1 -1
end if
ic21 = ic11 + 1
if (ps(igrid_working)%x(ic1,ic2,ic3,2) .lt. y(2)) then
ic12 = ic2
else
ic12 = ic2 -1
end if
ic22 = ic12 + 1
if (ps(igrid_working)%x(ic1,ic2,ic3,3) .lt. y(3)) then
ic13 = ic3
else
ic13 = ic3 -1
end if
ic23 = ic13 + 1
int_factor =0.5d0
if (ic11 .le. ixGlo1-2 .or. ic21 .ge. ixGhi1+2) then
int_factor = 0.05d0
end if
if (ic12 .le. ixGlo2-2 .or. ic22 .ge. ixGhi2+2) then
int_factor = 0.05d0
end if
if (ic13 .le. ixGlo3-2 .or. ic23 .ge. ixGhi3+2) then
int_factor = 0.05d0
end if
if (ic11 .eq. ixGlo1-1 .or. ic21 .eq. ixGhi1+1) then
int_factor = 0.1d0
end if
if (ic12 .eq. ixGlo2-1 .or. ic22 .eq. ixGhi2+1) then
int_factor = 0.1d0
end if
if (ic13 .eq. ixGlo3-1 .or. ic23 .eq. ixGhi3+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