24 integer,
parameter :: rk4=1, ark4=2
37 if (physics_type/=
'mhd')
call mpistop(
"GCA particles need magnetic field!")
38 if (
ndir/=3)
call mpistop(
"GCA particles need ndir=3!")
98 case(
'RK4',
'Rk4',
'rk4')
100 case(
'ARK4',
'ARk4',
'Ark4',
'ark4')
115 double precision :: b(
ndir), u(
ndir), magmom
116 double precision :: bnorm, lfac, vnorm, vperp, vpar
124 integer :: igrid, ipe_particle
125 integer :: n, idir, nparticles_local
133 rrd(n,idir) =
rng%unif_01()
137 {^
d&x(^
d,n) = xprobmin^
d + rrd(n+1,^
d) * (xprobmax^
d - xprobmin^
d)\}
159 if(ipe_particle ==
mype)
then
170 nparticles_local = nparticles_local + 1
186 vnorm = norm2(v(:, n))
187 vpar = sum(v(:, n) * b/bnorm)
188 vperp = sqrt(vnorm**2 - vpar**2)
196 magmom = m(n) * (vperp * lfac)**2 / (2.0d0 * bnorm)
204 call gca_update_payload(igrid,x(:,n),
particle(n)%self%u(:),q(n),m(n),defpayload,
ndefpayload,0.d0)
222 subroutine gca_fill_gridvars
227 double precision,
dimension(ixG^T,1:ndir) :: beta
228 double precision,
dimension(ixG^T,1:nw) :: w,wold
229 double precision :: current(ixg^t,7-2*
ndir:3)
230 double precision,
dimension(ixG^T,1:ndir) :: ve, bhat
231 double precision,
dimension(ixG^T) :: kappa, kappa_b, absb, tmp
232 integer :: igrid, iigrid, idir, idim
237 do iigrid=1,igridstail; igrid=igrids(iigrid);
238 w(ixg^t,1:nw) = ps(igrid)%w(ixg^t,1:nw)
239 call phys_to_primitive(ixg^
ll,ixg^
ll,w,ps(igrid)%x)
242 absb(ixg^t) = sqrt(sum(gridvars(igrid)%w(ixg^t,
bp(:))**2,dim=
ndim+1))
243 ve(ixg^t,1) = gridvars(igrid)%w(ixg^t,
ep(2)) * gridvars(igrid)%w(ixg^t,
bp(3)) &
244 - gridvars(igrid)%w(ixg^t,
ep(3)) * gridvars(igrid)%w(ixg^t,
bp(2))
245 ve(ixg^t,2) = gridvars(igrid)%w(ixg^t,
ep(3)) * gridvars(igrid)%w(ixg^t,
bp(1)) &
246 - gridvars(igrid)%w(ixg^t,
ep(1)) * gridvars(igrid)%w(ixg^t,
bp(3))
247 ve(ixg^t,3) = gridvars(igrid)%w(ixg^t,
ep(1)) * gridvars(igrid)%w(ixg^t,
bp(2)) &
248 - gridvars(igrid)%w(ixg^t,
ep(2)) * gridvars(igrid)%w(ixg^t,
bp(1))
250 where (absb(ixg^t) .gt. 0.d0)
251 ve(ixg^t,idir) = ve(ixg^t,idir) / absb(ixg^t)**2
253 ve(ixg^t,idir) = 0.d0
257 call mpistop(
"GCA FILL GRIDVARS: vE>c! ABORTING...")
259 if (any(ve .ne. ve))
then
260 call mpistop(
"GCA FILL GRIDVARS: NaNs IN vE! ABORTING...")
264 kappa(ixg^t) = 1.d0/sqrt(1.0d0 - sum(ve(ixg^t,:)**2,dim=
ndim+1)/
c_norm**2)
268 kappa_b(ixg^t) = absb(ixg^t) / kappa(ixg^t)
270 if (any(kappa_b .ne. kappa_b))
then
271 call mpistop(
"GCA FILL GRIDVARS: NaNs IN kappa_B! ABORTING...")
277 gridvars(igrid)%w(ixg^t,
grad_kappa_b(idim)) = tmp(ixg^t)
282 where (absb(ixg^t) .gt. 0.d0)
283 bhat(ixg^t,idir) = gridvars(igrid)%w(ixg^t,
bp(idir)) / absb(ixg^t)
285 bhat(ixg^t,idir) = 0.d0
288 if (any(bhat .ne. bhat))
then
289 call mpistop(
"GCA FILL GRIDVARS: NaNs IN bhat! ABORTING...")
295 call gradient(bhat(ixg^t,idir),ixg^
ll,ixg^
ll^lsub1,idim,tmp)
297 + bhat(ixg^t,idim) * tmp(ixg^t)
299 + ve(ixg^t,idim) * tmp(ixg^t)
300 call gradient(ve(ixg^t,idir),ixg^
ll,ixg^
ll^lsub1,idim,tmp)
302 + bhat(ixg^t,idim) * tmp(ixg^t)
304 + ve(ixg^t,idim) * tmp(ixg^t)
310 end subroutine gca_fill_gridvars
312 subroutine gca_integrate_particles(end_time)
316 double precision,
intent(in) :: end_time
318 double precision :: lfac, abss
321 double precision :: dt_p, tloc, y(
ndir+2),dydt(
ndir+2),ytmp(
ndir+2), euler_cfl, int_factor
323 double precision,
dimension(1:ndir) :: x, ve, e, b, bhat, x_new, vfluid, current
324 double precision,
dimension(1:ndir) :: drift1, drift2
325 double precision,
dimension(1:ndir) :: drift3, drift4, drift5, drift6, drift7
326 double precision,
dimension(1:ndir) :: bdotgradb, vedotgradb, gradkappab
327 double precision,
dimension(1:ndir) :: bdotgradve, vedotgradve
328 double precision,
dimension(1:ndir) :: gradbdrift, reldrift, bdotgradbdrift
329 double precision,
dimension(1:ndir) :: vedotgradbdrift, bdotgradvedrift
330 double precision,
dimension(1:ndir) :: vedotgradvedrift
331 double precision :: kappa, mr, upar, m, absb, gamma, q, mompar, vpar, veabs
332 double precision :: gradbdrift_abs, reldrift_abs, epar
333 double precision :: bdotgradbdrift_abs, vedotgradbdrift_abs
334 double precision :: bdotgradvedrift_abs, vedotgradvedrift_abs
335 double precision :: momentumpar1, momentumpar2, momentumpar3, momentumpar4
337 double precision,
parameter :: eps=1.0
d-6
339 double precision :: h1, hmin, h_old
340 integer :: nok, nbad, ic1^
d, ic2^
d, ierror, nvar
341 integer :: ipart, iipart, seed, ic^
d,igrid_particle, ipe_particle, ipe_working
342 logical :: int_choice
343 logical :: bc_applied
350 dt_p = gca_get_particle_dt(
particle(ipart), end_time)
365 call derivs_gca_rk(tk,y,k1)
366 tk = tloc + dt_p/2.d0
367 y(1:
ndir) = x + dt_p*k1(1:
ndir)/2.d0
369 call derivs_gca_rk(tk,y,k2)
370 tk = tloc + dt_p/2.d0
371 y(1:
ndir) = x + dt_p*k2(1:
ndir)/2.d0
373 call derivs_gca_rk(tk,y,k3)
377 call derivs_gca_rk(tk,y,k4)
380 y = y + dt_p/6.d0*(k1 + 2.d0*k2 + 2.d0*k3 + k4)
396 call derivs_gca(
particle(ipart)%self%time,y,dydt)
404 y(1:
ndir+2) = y(1:
ndir+2) + euler_cfl * dt_p * dydt(1:
ndir+2)
459 h1 = dt_p/2.0d0; hmin=1.0
d-9; h_old=dt_p/2.0d0
466 if (any(y .ne. y))
then
467 call mpistop(
"NaNs DETECTED IN GCA_INTEGRATE BEFORE ODEINT CALL! ABORTING...")
471 call odeint(y,nvar,tloc,tloc+dt_p,eps,h1,hmin,nok,nbad,derivs_gca_rk,
rkqs,ierror)
473 if (ierror /= 0)
then
474 print *,
"odeint returned error code", ierror
475 print *,
"1 means hmin too small, 2 means MAXSTP exceeded"
476 print *,
"Having a problem with particle", iipart
479 if (any(y .ne. y))
then
480 call mpistop(
"NaNs DETECTED IN GCA_INTEGRATE AFTER ODEINT CALL! ABORTING...")
503 absb = sqrt(sum(b(:)**2))
506 epar = sum(e(:)*bhat(:))
508 call cross(e,bhat,ve)
511 veabs = sqrt(sum(ve(:)**2))
513 kappa = 1.d0/sqrt(1.0d0 - sum(ve(:)**2)/
c_norm**2)
519 gamma = sqrt(1.0d0+upar**2/
c_norm**2+2.0d0*mr*absb/m/
c_norm**2)*kappa
541 end subroutine gca_integrate_particles
543 subroutine derivs_gca_rk(t_s,y,dydt)
546 double precision :: t_s, y(
ndir+2)
547 double precision :: dydt(
ndir+2)
549 double precision,
dimension(ndir):: ve, b, e, x, bhat, bdotgradb, vedotgradb, gradkappab, vfluid, current
550 double precision,
dimension(ndir):: bdotgradve, vedotgradve, u, utmp1, utmp2, utmp3
551 double precision :: upar, mr, gamma, absb, q, m, epar, kappa
566 if (any(x .ne. x))
then
567 write(*,*)
"ERROR IN DERIVS_GCA_RK: NaNs IN X OR Y!"
570 call mpistop(
"ABORTING...")
587 if (any(b .ne. b) .or. any(e .ne. e) &
588 .or. any(bdotgradb .ne. bdotgradb) .or. any(vedotgradb .ne. vedotgradb) &
589 .or. any(gradkappab .ne. gradkappab) .or. any(bdotgradve .ne. bdotgradve) &
590 .or. any(vedotgradve .ne. vedotgradve))
then
591 write(*,*)
"ERROR IN DERIVS_GCA_RK: NaNs IN FIELD QUANTITIES!"
594 write(*,*)
"bdotgradb",bdotgradb
595 write(*,*)
"vEdotgradb",vedotgradb
596 write(*,*)
"gradkappaB",gradkappab
597 write(*,*)
"bdotgradvE",bdotgradve
598 write(*,*)
"vEdotgradvE",vedotgradve
599 call mpistop(
"ABORTING...")
602 absb = sqrt(sum(b(:)**2))
603 if (absb .gt. 0.d0)
then
608 epar = sum(e(:)*bhat(:))
610 call cross(e,bhat,ve)
611 if (absb .gt. 0.d0)
then
618 kappa = 1.d0/sqrt(1.0d0 - sum(ve(:)**2)/
c_norm**2)
619 gamma = sqrt(1.0d0+upar**2/
c_norm**2+2.0d0*mr*absb/m/
c_norm**2)*kappa
625 if (absb .gt. 0.d0)
then
626 utmp1(1:
ndir) = bhat(1:
ndir)/(absb/kappa**2)
630 utmp2(1:
ndir) = mr/(gamma*q)*gradkappab(1:
ndir) &
631 + m/q* (upar**2/gamma*bdotgradb(1:
ndir) + upar*vedotgradb(1:
ndir) &
632 + upar*bdotgradve(1:
ndir) + gamma*vedotgradve(1:
ndir))
634 utmp2(1:
ndir) = utmp2(1:
ndir) + upar*epar/(gamma)*ve(1:
ndir)
637 call cross(utmp1,utmp2,utmp3)
641 dydt(1:
ndir) = ( u(1:
ndir) + upar/gamma * bhat(1:
ndir) )
642 dydt(
ndir+1) = q/m*epar - mr/(m*gamma) * sum(bhat(:)*gradkappab(:)) &
643 + sum(ve(:)*(upar*bdotgradb(:)+gamma*vedotgradb(:)))
646 end subroutine derivs_gca_rk
648 subroutine derivs_gca(t_s,y,dydt)
651 double precision :: t_s, y(
ndir+2)
652 double precision :: dydt(
ndir+2)
654 double precision,
dimension(ndir):: ve, b, e, x, bhat, bdotgradb, vedotgradb, gradkappab, vfluid, current
655 double precision,
dimension(ndir):: bdotgradve, vedotgradve, u, utmp1, utmp2, utmp3
656 double precision :: upar, mr, gamma, absb, q, m, epar, kappa
668 if (any(x .ne. x))
then
669 call mpistop(
"ERROR IN DERIVS_GCA: NaNs IN X! ABORTING...")
686 absb = sqrt(sum(b(:)**2))
687 if (absb .gt. 0.d0)
then
693 epar = sum(e(:)*bhat(:))
694 call cross(e,bhat,ve)
695 if (absb .gt. 0.d0) ve(1:
ndir) = ve(1:
ndir) / absb
698 kappa = 1.d0/sqrt(1.0d0 - sum(ve(:)**2)/
c_norm**2)
699 gamma = sqrt(1.0d0+upar**2/
c_norm**2+2.0d0*mr*absb/m/
c_norm**2)*kappa
704 if (absb .gt. 0.d0)
then
705 utmp1(1:
ndir) = bhat(1:
ndir)/(absb/kappa**2)
709 utmp2(1:
ndir) = mr/(gamma*q)*gradkappab(1:
ndir) &
710 + m/q* (upar**2/gamma*bdotgradb(1:
ndir) + upar*vedotgradb(1:
ndir) &
711 + upar*bdotgradve(1:
ndir) + gamma*vedotgradve(1:
ndir))
713 utmp2(1:
ndir) = utmp2(1:
ndir) + upar*epar/(gamma)*ve(1:
ndir)
716 call cross(utmp1,utmp2,utmp3)
720 dydt(1:
ndir) = ( u(1:
ndir) + upar/gamma * bhat(1:
ndir) )
721 dydt(
ndir+1) = q/m*epar - mr/(m*gamma) * sum(bhat(:)*gradkappab(:)) &
722 + sum(ve(:)*(upar*bdotgradb(:)+gamma*vedotgradb(:)))
725 end subroutine derivs_gca
728 subroutine gca_update_payload(igrid,xpart,upart,qpart,mpart,mypayload,mynpayload,particle_time)
730 integer,
intent(in) :: igrid,mynpayload
731 double precision,
intent(in) :: xpart(1:
ndir),upart(1:
ndir),qpart,mpart,particle_time
732 double precision,
intent(out) :: mypayload(mynpayload)
733 double precision,
dimension(1:ndir) :: ve, e, b, bhat, vfluid, current
734 double precision,
dimension(1:ndir) :: drift1, drift2
735 double precision,
dimension(1:ndir) :: drift3, drift4, drift5, drift6, drift7
736 double precision,
dimension(1:ndir) :: bdotgradb, vedotgradb, gradkappab
737 double precision,
dimension(1:ndir) :: bdotgradve, vedotgradve
738 double precision,
dimension(1:ndir) :: gradbdrift, reldrift, bdotgradbdrift
739 double precision,
dimension(1:ndir) :: vedotgradbdrift, bdotgradvedrift
740 double precision,
dimension(1:ndir) :: vedotgradvedrift
741 double precision :: kappa, upar, absb, gamma, vpar, veabs
742 double precision :: gradbdrift_abs, reldrift_abs, epar
743 double precision :: bdotgradbdrift_abs, vedotgradbdrift_abs
744 double precision :: bdotgradvedrift_abs, vedotgradvedrift_abs
745 double precision :: momentumpar1, momentumpar2, momentumpar3, momentumpar4
756 absb = sqrt(sum(b(:)**2))
758 epar = sum(e(:)*bhat(:))
759 call cross(e,bhat,ve)
761 veabs = sqrt(sum(ve(:)**2))
763 kappa = 1.d0/sqrt(1.0d0 - sum(ve(:)**2)/
c_norm**2)
767 vpar = upart(1)/upart(3)
776 drift1(1:
ndir) = bhat(1:
ndir)/(absb/kappa**2)
777 drift2(1:
ndir) = upart(2)/(upart(3)*qpart)*gradkappab(1:
ndir)
779 call cross(drift1,drift2,gradbdrift)
780 gradbdrift_abs = sqrt(sum(gradbdrift(:)**2))
782 drift3(1:
ndir) = upar*epar/upart(3)*ve(1:
ndir)
783 call cross(drift1,drift3,reldrift)
784 reldrift_abs = sqrt(sum(reldrift(:)**2))
786 drift4(1:
ndir) = mpart/qpart* ( upar**2/upart(3)*bdotgradb(1:
ndir))
787 call cross(drift1,drift4,bdotgradbdrift)
788 bdotgradbdrift_abs = sqrt(sum(bdotgradbdrift(:)**2))
790 drift5(1:
ndir) = mpart/qpart* ( upar*vedotgradb(1:
ndir))
791 call cross(drift1,drift5,vedotgradbdrift)
792 vedotgradbdrift_abs = sqrt(sum(vedotgradbdrift(:)**2))
794 drift6(1:
ndir) = mpart/qpart* ( upar*bdotgradve(1:
ndir))
795 call cross(drift1,drift6,bdotgradvedrift)
796 bdotgradvedrift_abs = sqrt(sum(bdotgradvedrift(:)**2))
798 drift7(1:
ndir) = mpart/qpart* (upart(3)*vedotgradve(1:
ndir))
799 call cross(drift1,drift7,vedotgradvedrift)
800 vedotgradvedrift_abs = sqrt(sum(vedotgradvedrift(:)**2))
802 momentumpar1 = qpart/mpart*epar
803 momentumpar2 = -(upart(2)/mpart/upart(3))*sum(bhat(:)*gradkappab(:))
804 momentumpar3 = upar*sum(ve(:)*bdotgradb(:))
805 momentumpar4 = upart(3)*sum(ve(:)*vedotgradb(:))
808 if (mynpayload > 0)
then
810 mypayload(1) = sqrt(2.0d0*mpart*upart(2)*absb)/abs(qpart*absb)
812 if (mynpayload > 1)
then
814 mypayload(2) = atan2(sqrt((2.0d0*upart(2)*absb)/(mpart*upart(3)**2)),vpar)
816 if (mynpayload > 2)
then
818 mypayload(3) = sqrt((2.0d0*upart(2)*absb)/(mpart*upart(3)**2))
820 if (mynpayload > 3)
then
822 mypayload(4) = momentumpar1
824 if (mynpayload > 4)
then
826 mypayload(5) = momentumpar2
828 if (mynpayload > 5)
then
830 mypayload(6) = momentumpar3
832 if (mynpayload > 6)
then
834 mypayload(7) = momentumpar4
836 if (mynpayload > 7)
then
840 if (mynpayload > 8)
then
842 mypayload(9) = gradbdrift_abs
844 if (mynpayload > 9)
then
846 mypayload(10) = reldrift_abs
848 if (mynpayload > 10)
then
850 mypayload(11) = bdotgradbdrift_abs
852 if (mynpayload > 11)
then
854 mypayload(12) = vedotgradbdrift_abs
856 if (mynpayload > 12)
then
858 mypayload(13) = bdotgradvedrift_abs
860 if (mynpayload > 13)
then
862 mypayload(14) = vedotgradvedrift_abs
865 end subroutine gca_update_payload
867 function gca_get_particle_dt(partp, end_time)
result(dt_p)
871 double precision,
intent(in) :: end_time
872 double precision :: dt_p
874 double precision :: tout, dt_particles_mype, dt_cfl0, dt_cfl1, dt_a
875 double precision :: dxmin,
vp, a, gammap
877 double precision :: ap0, ap1, dt_cfl_ap0, dt_cfl_ap1, dt_cfl_ap
878 double precision :: dt_euler, dt_tmp
880 double precision :: cfl, uparcfl
881 double precision,
parameter :: uparmin=1.0
d-6*const_c
882 integer :: ipart, iipart, nout, ic^
d, igrid_particle, ipe_particle, ipe
883 logical :: bc_applied
895 dt_tmp = (end_time - partp%self%time)
896 if(dt_tmp .le. 0.0d0) dt_tmp = smalldouble
904 y(
ndir+1) = partp%self%u(1)
905 y(
ndir+2) = partp%self%u(2)
909 call derivs_gca(partp%self%time,y,dydt)
915 vp = sqrt(sum(v(:)**2))
917 dt_cfl0 = dxmin / max(
vp, smalldouble)
918 dt_cfl_ap0 = uparcfl * abs(max(abs(y(
ndir+1)),uparmin) / max(abs(ap0), smalldouble))
923 dt_euler = min(dt_tmp,dt_cfl0,dt_cfl_ap0)
927 partp%self%u(1) = y(
ndir+1)
928 partp%self%u(2) = y(
ndir+2)
935 call derivs_gca_rk(partp%self%time+dt_euler,y,dydt)
943 vp = sqrt(sum(v(:)**2))
945 dt_cfl1 = dxmin / max(
vp, smalldouble)
946 dt_cfl_ap1 = uparcfl * abs(max(abs(y(
ndir+1)),uparmin) / max(abs(ap1), smalldouble))
949 dt_tmp = min(dt_euler, dt_cfl1, dt_cfl_ap1)
952 partp%self%u(1) = ytmp(
ndir+1)
953 partp%self%u(2) = ytmp(
ndir+2)
976 end function gca_get_particle_dt
Module with geometry-related routines (e.g., divergence, curl)
subroutine gradient(q, ixil, ixol, idir, gradq, nth_in)
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer ixghi
Upper index of grid block arrays.
double precision global_time
The global simulation time.
integer, parameter ndim
Number of spatial dimensions for grid variables.
integer, parameter rpxmin
integer icomm
The MPI communicator.
integer mype
The rank of the current MPI task.
double precision, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ierrmpi
A global MPI error return code.
double precision c_norm
Normalised speed of light.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter ixglo
Lower index of grid block arrays (always 1)
This module packages odeint from numerical recipes.
subroutine odeint(ystart, nvar, x1, x2, eps, h1, hmin, nok, nbad, derivs, rkqs, ierror)
subroutine rkqs(y, dydx, n, x, htry, eps, yscal, hdid, hnext, derivs)
Module with shared functionality for all the particle movers.
pure subroutine limit_dt_endtime(t_left, dt_p)
integer num_particles
Number of particles.
double precision const_dt_particles
If positive, a constant time step for the particles.
subroutine write_particle_output()
logical function particle_in_igrid(ipart, igrid)
Quick check if particle is still in igrid.
double precision dtsave_particles
Time interval to save particles.
integer ngridvars
Number of variables for grid field.
subroutine get_efield(igrid, x, tloc, e)
integer nusrpayload
Number of user-defined payload variables for a particle.
subroutine find_particle_ipe(x, igrid_particle, ipe_particle)
pure subroutine get_lfac_from_velocity(v, lfac)
Get Lorentz factor from velocity.
integer integrator
Integrator to be used for particles.
integer, dimension(:), allocatable ep
Variable index for electric field.
subroutine fill_gridvars_default
This routine fills the particle grid variables with the default for mhd, i.e. only E and B.
subroutine get_bfield(igrid, x, tloc, b)
integer npayload
Number of total payload variables for a particle.
integer, dimension(:), allocatable particles_active_on_mype
Array of identity numbers of active particles in current processor.
logical relativistic
Use a relativistic particle mover?
character(len=name_len) integrator_type_particles
String describing the particle integrator type.
subroutine get_vec(ix, igrid, x, tloc, vec)
subroutine push_particle_into_particles_on_mype(ipart)
integer nparticles_active_on_mype
Number of active particles in current processor.
procedure(sub_define_additional_gridvars), pointer particles_define_additional_gridvars
integer nparticles
Identity number and total number of particles.
integer, dimension(:), allocatable bp
Variable index for magnetic field.
type(particle_ptr), dimension(:), allocatable particle
integer, dimension(:), allocatable vp
Variable index for fluid velocity.
procedure(sub_integrate), pointer particles_integrate
subroutine cross(a, b, c)
procedure(sub_fill_gridvars), pointer particles_fill_gridvars
double precision t_next_output
Time to write next particle output.
integer ndefpayload
Number of default payload variables for a particle.
integer igrid_working
set the current igrid for the particle integrator:
integer, dimension(:), allocatable jp
Variable index for current.
integer ipart_working
set the current ipart for the particle integrator:
double precision particles_cfl
Particle CFL safety factor.
integer rhop
Variable index for density.
Particle mover with Newtonian/relativistic Guiding Center Approximation (GCA) By Jannis Teunissen,...
integer, dimension(:), allocatable, public, protected b_dot_grad_b
Variable index for (B . grad)B (curvature B drift)
subroutine, public gca_init()
subroutine, public gca_create_particles()
integer, dimension(:), allocatable, public, protected b_dot_grad_ve
Variable index for polarization drift.
integer, dimension(:), allocatable, public, protected ve_dot_grad_ve
Variable index for polarization drift.
integer, dimension(:), allocatable, public, protected grad_kappa_b
Variable index for gradient B, with relativistic correction 1/kappa where kappa = 1/sqrt(1 - E_perp^2...
integer, dimension(:), allocatable, public, protected ve_dot_grad_b
Variable index for curvature drift.
Module with all the methods that users can customize in AMRVAC.
procedure(check_particle), pointer usr_check_particle
procedure(create_particles), pointer usr_create_particles
procedure(update_payload), pointer usr_update_payload
procedure(particle_fields), pointer usr_particle_fields