MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_particle_gca.t
Go to the documentation of this file.
1 !> Particle mover with Newtonian/relativistic Guiding Center Approximation (GCA)
2 !> By Jannis Teunissen, Bart Ripperda, Oliver Porth, and Fabio Bacchini (2016-2020)
5 
6  private
7 
8  !> Variable index for gradient B, with relativistic correction 1/kappa
9  !> where kappa = 1/sqrt(1 - E_perp^2/B^2)
10  integer, protected, allocatable :: grad_kappa_b(:)
11  !> Variable index for (B . grad)B (curvature B drift)
12  integer, protected, allocatable :: b_dot_grad_b(:)
13 
14  ! ExB related drifts (vE = ExB/B^2)
15  !> Variable index for curvature drift
16  integer, protected, allocatable :: ve_dot_grad_b(:)
17  !> Variable index for polarization drift
18  integer, protected, allocatable :: b_dot_grad_ve(:)
19  !> Variable index for polarization drift
20  integer, protected, allocatable :: ve_dot_grad_ve(:)
21 
22  public :: gca_init
23  public :: gca_create_particles
24  integer, parameter :: rk4=1, ark4=2
25 
26  ! Variables
27  public :: bp, ep, grad_kappa_b, b_dot_grad_b
29  public :: vp, jp, rhop
30 
31 contains
32 
33  subroutine gca_init()
35  integer :: idir, nwx
36 
37  if (physics_type/='mhd') call mpistop("GCA particles need magnetic field!")
38  if (ndir/=3) call mpistop("GCA particles need ndir=3!")
39 
40  nwx = 0
41  allocate(bp(ndir))
42  do idir = 1, ndir
43  nwx = nwx + 1
44  bp(idir) = nwx
45  end do
46  allocate(ep(ndir))
47  do idir = 1, ndir
48  nwx = nwx + 1
49  ep(idir) = nwx
50  end do
51  allocate(grad_kappa_b(ndir))
52  do idir = 1, ndir
53  nwx = nwx + 1
54  grad_kappa_b(idir) = nwx
55  end do
56  allocate(b_dot_grad_b(ndir))
57  do idir = 1, ndir
58  nwx = nwx + 1
59  b_dot_grad_b(idir) = nwx
60  end do
61  allocate(ve_dot_grad_b(ndir))
62  do idir = 1, ndir
63  nwx = nwx + 1
64  ve_dot_grad_b(idir) = nwx
65  end do
66  allocate(b_dot_grad_ve(ndir))
67  do idir = 1, ndir
68  nwx = nwx + 1
69  b_dot_grad_ve(idir) = nwx
70  end do
71  allocate(ve_dot_grad_ve(ndir))
72  do idir = 1, ndir
73  nwx = nwx + 1
74  ve_dot_grad_ve(idir) = nwx
75  end do
76  allocate(vp(ndir))
77  do idir = 1, ndir
78  nwx = nwx + 1
79  vp(idir) = nwx
80  end do
81  allocate(jp(ndir))
82  do idir = 1, ndir
83  nwx = nwx + 1
84  jp(idir) = nwx
85  end do
86  nwx = nwx + 1 ! density
87  rhop = nwx
88 
89  ngridvars=nwx
90 
92 
93  if (associated(particles_define_additional_gridvars)) then
95  end if
96 
97  select case(integrator_type_particles)
98  case('RK4','Rk4','rk4')
99  integrator = rk4
100  case('ARK4','ARk4','Ark4','ark4')
101  integrator = ark4
102  case default
103  integrator = ark4
104  end select
105 
107 
108  end subroutine gca_init
109 
110  subroutine gca_create_particles()
111  ! initialise the particles
114 
115  double precision :: b(ndir), u(ndir), magmom
116  double precision :: bnorm, lfac, vnorm, vperp, vpar
117  integer :: igrid, ipe_particle
118  integer :: n, idir, nparticles_local
119  double precision :: x(3, num_particles)
120  double precision :: v(3, num_particles)
121  double precision :: q(num_particles)
122  double precision :: m(num_particles)
123  double precision :: rrd(num_particles,ndir)
124  double precision :: defpayload(ndefpayload)
125  double precision :: usrpayload(nusrpayload)
126  logical :: follow(num_particles), check
127 
128  if (mype==0) then
129  if (.not. associated(usr_create_particles)) then
130  ! Randomly distributed
131  do idir=1,ndir
132  do n = 1, num_particles
133  rrd(n,idir) = rng%unif_01()
134  end do
135  end do
136  do n=1, num_particles
137  {^d&x(^d,n) = xprobmin^d + rrd(n+1,^d) * (xprobmax^d - xprobmin^d)\}
138  end do
139  else
140  call usr_create_particles(num_particles, x, v, q, m, follow)
141  end if
142  end if
143 
144  call mpi_bcast(x,3*num_particles,mpi_double_precision,0,icomm,ierrmpi)
145  call mpi_bcast(v,3*num_particles,mpi_double_precision,0,icomm,ierrmpi)
146  call mpi_bcast(q,num_particles,mpi_double_precision,0,icomm,ierrmpi)
147  call mpi_bcast(m,num_particles,mpi_double_precision,0,icomm,ierrmpi)
148  call mpi_bcast(follow,num_particles,mpi_logical,0,icomm,ierrmpi)
149 
150  nparticles_local = 0
151 
152  ! first find ipe and igrid responsible for particle
153  do n = 1, num_particles
154  call find_particle_ipe(x(:, n),igrid,ipe_particle)
155 
156  particle(n)%igrid = igrid
157  particle(n)%ipe = ipe_particle
158 
159  if(ipe_particle == mype) then
160  check = .true.
161 
162  ! Check for user-defined modifications or rejection conditions
163  if (associated(usr_check_particle)) call usr_check_particle(igrid, x(:,n), v(:,n), q(n), m(n), follow(n), check)
164  if (check) then
166  else
167  cycle
168  end if
169 
170  nparticles_local = nparticles_local + 1
171 
172  call get_lfac_from_velocity(v(:, n), lfac)
173 
174  allocate(particle(n)%self)
175  particle(n)%self%x = x(:, n)
176  particle(n)%self%q = q(n)
177  particle(n)%self%m = m(n)
178  particle(n)%self%follow = follow(n)
179  particle(n)%self%index = n
180  particle(n)%self%time = global_time
181  particle(n)%self%dt = 0.0d0
182 
183  call get_vec(bp, igrid, x(:, n), particle(n)%self%time, b)
184 
185  bnorm = norm2(b(:))
186  vnorm = norm2(v(:, n))
187  vpar = sum(v(:, n) * b/bnorm)
188  vperp = sqrt(vnorm**2 - vpar**2)
189 
190  ! The momentum vector u(1:3) is filled with the following components
191 
192  ! parallel momentum component (gamma v||)
193  particle(n)%self%u(1) = lfac * vpar
194 
195  ! Mr: the conserved magnetic moment
196  magmom = m(n) * (vperp * lfac)**2 / (2.0d0 * bnorm)
197  particle(n)%self%u(2) = magmom
198 
199  ! Lorentz factor
200  particle(n)%self%u(3) = lfac
201 
202  ! initialise payloads for GCA module
203  allocate(particle(n)%payload(npayload))
204  call gca_update_payload(igrid,x(:,n),particle(n)%self%u(:),q(n),m(n),defpayload,ndefpayload,0.d0)
205  particle(n)%payload(1:ndefpayload) = defpayload
206  if (associated(usr_update_payload)) then
207  call usr_update_payload(igrid,x(:,n),particle(n)%self%u(:),q(n),m(n),usrpayload,nusrpayload,0.d0)
208  particle(n)%payload(ndefpayload+1:npayload) = usrpayload
209  end if
210  end if
211  end do
212 
213  call mpi_allreduce(nparticles_local,nparticles,1,mpi_integer,mpi_sum,icomm,ierrmpi)
214 
215  end subroutine gca_create_particles
216 
217  subroutine gca_fill_gridvars
220  use mod_geometry
221 
222  integer :: igrid, iigrid, idir, idim
223  double precision, dimension(ixG^T,1:ndir) :: beta
224  double precision, dimension(ixG^T,1:nw) :: w,wold
225  double precision :: current(ixG^T,7-2*ndir:3)
226  integer :: idirmin
227  double precision, dimension(ixG^T,1:ndir) :: vE, bhat
228  double precision, dimension(ixG^T) :: kappa, kappa_B, absB, tmp
229 
230  call fill_gridvars_default()
231 
232  do iigrid=1,igridstail; igrid=igrids(iigrid);
233  w(ixg^t,1:nw) = ps(igrid)%w(ixg^t,1:nw)
234  call phys_to_primitive(ixg^ll,ixg^ll,w,ps(igrid)%x)
235 
236  ! grad(kappa B)
237  absb(ixg^t) = sqrt(sum(gridvars(igrid)%w(ixg^t,bp(:))**2,dim=ndim+1))
238  ve(ixg^t,1) = gridvars(igrid)%w(ixg^t,ep(2)) * gridvars(igrid)%w(ixg^t,bp(3)) &
239  - gridvars(igrid)%w(ixg^t,ep(3)) * gridvars(igrid)%w(ixg^t,bp(2))
240  ve(ixg^t,2) = gridvars(igrid)%w(ixg^t,ep(3)) * gridvars(igrid)%w(ixg^t,bp(1)) &
241  - gridvars(igrid)%w(ixg^t,ep(1)) * gridvars(igrid)%w(ixg^t,bp(3))
242  ve(ixg^t,3) = gridvars(igrid)%w(ixg^t,ep(1)) * gridvars(igrid)%w(ixg^t,bp(2)) &
243  - gridvars(igrid)%w(ixg^t,ep(2)) * gridvars(igrid)%w(ixg^t,bp(1))
244  do idir=1,ndir
245  where (absb(ixg^t) .gt. 0.d0)
246  ve(ixg^t,idir) = ve(ixg^t,idir) / absb(ixg^t)**2
247  elsewhere
248  ve(ixg^t,idir) = 0.d0
249  end where
250  end do
251  if (any(sum(ve(ixg^t,:)**2,dim=ndim+1) .ge. c_norm**2) .and. relativistic) then
252  call mpistop("GCA FILL GRIDVARS: vE>c! ABORTING...")
253  end if
254  if (any(ve .ne. ve)) then
255  call mpistop("GCA FILL GRIDVARS: NaNs IN vE! ABORTING...")
256  end if
257 
258  if (relativistic) then
259  kappa(ixg^t) = 1.d0/sqrt(1.0d0 - sum(ve(ixg^t,:)**2,dim=ndim+1)/c_norm**2)
260  else
261  kappa(ixg^t) = 1.d0
262  end if
263  kappa_b(ixg^t) = absb(ixg^t) / kappa(ixg^t)
264 
265  if (any(kappa_b .ne. kappa_b)) then
266  call mpistop("GCA FILL GRIDVARS: NaNs IN kappa_B! ABORTING...")
267  end if
268 
269  tmp=0.d0
270  do idim=1,ndim
271  call gradient(kappa_b,ixg^ll,ixg^ll^lsub1,idim,tmp)
272  gridvars(igrid)%w(ixg^t,grad_kappa_b(idim)) = tmp(ixg^t)
273  end do
274 
275  ! bhat
276  do idir=1,ndir
277  where (absb(ixg^t) .gt. 0.d0)
278  bhat(ixg^t,idir) = gridvars(igrid)%w(ixg^t,bp(idir)) / absb(ixg^t)
279  elsewhere
280  bhat(ixg^t,idir) = 0.d0
281  end where
282  end do
283  if (any(bhat .ne. bhat)) then
284  call mpistop("GCA FILL GRIDVARS: NaNs IN bhat! ABORTING...")
285  end if
286 
287  do idir=1,ndir
288  ! (b dot grad) b and the other directional derivatives
289  do idim=1,ndim
290  call gradient(bhat(ixg^t,idir),ixg^ll,ixg^ll^lsub1,idim,tmp)
291  gridvars(igrid)%w(ixg^t,b_dot_grad_b(idir)) = gridvars(igrid)%w(ixg^t,b_dot_grad_b(idir)) &
292  + bhat(ixg^t,idim) * tmp(ixg^t)
293  gridvars(igrid)%w(ixg^t,ve_dot_grad_b(idir)) = gridvars(igrid)%w(ixg^t,ve_dot_grad_b(idir)) &
294  + ve(ixg^t,idim) * tmp(ixg^t)
295  call gradient(ve(ixg^t,idir),ixg^ll,ixg^ll^lsub1,idim,tmp)
296  gridvars(igrid)%w(ixg^t,b_dot_grad_ve(idir)) = gridvars(igrid)%w(ixg^t,b_dot_grad_ve(idir)) &
297  + bhat(ixg^t,idim) * tmp(ixg^t)
298  gridvars(igrid)%w(ixg^t,ve_dot_grad_ve(idir)) = gridvars(igrid)%w(ixg^t,ve_dot_grad_ve(idir)) &
299  + ve(ixg^t,idim) * tmp(ixg^t)
300  end do
301  end do
302 
303  end do
304 
305  end subroutine gca_fill_gridvars
306 
307  subroutine gca_integrate_particles(end_time)
308  use mod_odeint
311  double precision, intent(in) :: end_time
312 
313  double precision :: lfac, absS
314  double precision :: defpayload(ndefpayload)
315  double precision :: usrpayload(nusrpayload)
316  double precision :: dt_p, tloc, y(ndir+2),dydt(ndir+2),ytmp(ndir+2), euler_cfl, int_factor
317  double precision :: tk, k1(ndir+2),k2(ndir+2),k3(ndir+2),k4(ndir+2)
318  double precision, dimension(1:ndir) :: x, vE, e, b, bhat, x_new, vfluid, current
319  double precision, dimension(1:ndir) :: drift1, drift2
320  double precision, dimension(1:ndir) :: drift3, drift4, drift5, drift6, drift7
321  double precision, dimension(1:ndir) :: bdotgradb, vEdotgradb, gradkappaB
322  double precision, dimension(1:ndir) :: bdotgradvE, vEdotgradvE
323  double precision, dimension(1:ndir) :: gradBdrift, reldrift, bdotgradbdrift
324  double precision, dimension(1:ndir) :: vEdotgradbdrift, bdotgradvEdrift
325  double precision, dimension(1:ndir) :: vEdotgradvEdrift
326  double precision :: kappa, Mr, upar, m, absb, gamma, q, mompar, vpar, vEabs
327  double precision :: gradBdrift_abs, reldrift_abs, epar
328  double precision :: bdotgradbdrift_abs, vEdotgradbdrift_abs
329  double precision :: bdotgradvEdrift_abs, vEdotgradvEdrift_abs
330  double precision :: momentumpar1, momentumpar2, momentumpar3, momentumpar4
331  ! Precision of time-integration:
332  double precision,parameter :: eps=1.0d-6
333  ! for odeint:
334  double precision :: h1, hmin, h_old
335  integer :: nok, nbad, ic1^D, ic2^D, ierror, nvar
336  integer :: ipart, iipart, seed, ic^D,igrid_particle, ipe_particle, ipe_working
337  logical :: int_choice
338  logical :: BC_applied
339 
340  nvar=ndir+2
341 
342  do iipart=1,nparticles_active_on_mype
343  ipart = particles_active_on_mype(iipart)
344  int_choice = .false.
345  dt_p = gca_get_particle_dt(particle(ipart), end_time)
346  particle(ipart)%self%dt = dt_p
347 
348  igrid_working = particle(ipart)%igrid
349  ipart_working = particle(ipart)%self%index
350  tloc = particle(ipart)%self%time
351  x(1:ndir) = particle(ipart)%self%x(1:ndir)
352 
353  select case (integrator)
354  case (rk4)
355  ! Runge-Kutta order 4
356  tk = tloc
357  y(1:ndir) = x(1:ndir) ! position of guiding center
358  y(ndir+1) = particle(ipart)%self%u(1) ! parallel momentum component (gamma v||)
359  y(ndir+2) = particle(ipart)%self%u(2) ! conserved magnetic moment Mr
360  call derivs_gca_rk(tk,y,k1)
361  tk = tloc + dt_p/2.d0
362  y(1:ndir) = x + dt_p*k1(1:ndir)/2.d0
363  y(ndir+1) = particle(ipart)%self%u(1) + dt_p*k1(ndir+1)/2.d0
364  call derivs_gca_rk(tk,y,k2)
365  tk = tloc + dt_p/2.d0
366  y(1:ndir) = x + dt_p*k2(1:ndir)/2.d0
367  y(ndir+1) = particle(ipart)%self%u(1) + dt_p*k2(ndir+1)/2.d0
368  call derivs_gca_rk(tk,y,k3)
369  tk = tloc + dt_p
370  y(1:ndir) = x + dt_p*k3(1:ndir)
371  y(ndir+1) = particle(ipart)%self%u(1) + dt_p*k3(ndir+1)
372  call derivs_gca_rk(tk,y,k4)
373  y(1:ndir) = x(1:ndir) ! position of guiding center
374  y(ndir+1) = particle(ipart)%self%u(1) ! parallel momentum component (gamma v||)
375  y = y + dt_p/6.d0*(k1 + 2.d0*k2 + 2.d0*k3 + k4)
376  particle(ipart)%self%x(1:ndir) = y(1:ndir)
377  particle(ipart)%self%u(1) = y(ndir+1)
378 
379  case (ark4)
380  ! Adaptive stepwidth RK4:
381  ! initial solution vector:
382  y(1:ndir) = x(1:ndir) ! position of guiding center
383  y(ndir+1) = particle(ipart)%self%u(1) ! parallel momentum component (gamma v||)
384  y(ndir+2) = particle(ipart)%self%u(2) ! conserved magnetic moment Mr
385  ! y(ndir+3) = particle(ipart)%self%u(3) ! Lorentz factor of particle
386 
387  ! we temporarily save the solution vector, to replace the one from the euler
388  ! timestep after euler integration
389  ytmp=y
390 
391  call derivs_gca(particle(ipart)%self%time,y,dydt)
392 
393  ! make an Euler step with the proposed timestep:
394  ! factor to ensure we capture all particles near the internal ghost cells.
395  ! Can be adjusted during a run, after an interpolation error.
396  euler_cfl=2.5d0
397 
398  ! new solution vector:
399  y(1:ndir+2) = y(1:ndir+2) + euler_cfl * dt_p * dydt(1:ndir+2)
400  particle(ipart)%self%x(1:ndir) = y(1:ndir) ! position of guiding center
401  particle(ipart)%self%u(1) = y(ndir+1) ! parallel momentum component(gamma v||)
402  particle(ipart)%self%u(2) = y(ndir+2) ! conserved magnetic moment
403 
404  ! check if the particle is in the internal ghost cells
405  int_factor =1.0d0
406 
408  ! if particle is not in the grid an euler timestep is taken instead of a RK4
409  ! timestep. Then based on that we do an interpolation and check how much further
410  ! the timestep for the RK4 has to be restricted.
411  ! factor to make integration more accurate for particles near the internal
412  ! ghost cells. This factor can be changed during integration after an
413  ! interpolation error. But one should be careful with timesteps for i/o
414 
415  ! flat interpolation:
416  {ic^d = int((y(^d)-rnode(rpxmin^d_,igrid_working))/rnode(rpdx^d_,igrid_working)) + 1 + nghostcells\}
417 
418  ! linear interpolation:
419  {
420  if (ps(igrid_working)%x({ic^dd},^d) .lt. y(^d)) then
421  ic1^d = ic^d
422  else
423  ic1^d = ic^d -1
424  end if
425  ic2^d = ic1^d + 1
426  \}
427 
428  int_factor =0.5d0
429 
430  {^d&
431  if (ic1^d .le. ixglo^d-2 .or. ic2^d .ge. ixghi^d+2) then
432  int_factor = 0.05d0
433  end if
434  \}
435 
436  {^d&
437  if (ic1^d .eq. ixglo^d-1 .or. ic2^d .eq. ixghi^d+1) then
438  int_factor = 0.1d0
439  end if
440  \}
441 
442  dt_p=int_factor*dt_p
443  end if
444 
445  ! replace the solution vector with the original as it was before the Euler timestep
446  y(1:ndir+2) = ytmp(1:ndir+2)
447 
448  particle(ipart)%self%x(1:ndir) = ytmp(1:ndir) ! position of guiding center
449  particle(ipart)%self%u(1) = ytmp(ndir+1) ! parallel momentum component (gamma v||)
450  particle(ipart)%self%u(2) = ytmp(ndir+2) ! conserved magnetic moment
451 
452  ! specify a minimum step hmin. If the timestep reaches this minimum, multiply by
453  ! a factor 100 to make sure the RK integration doesn't crash
454  h1 = dt_p/2.0d0; hmin=1.0d-9; h_old=dt_p/2.0d0
455 
456  if(h1 .lt. hmin)then
457  h1=hmin
458  dt_p=2.0d0*h1
459  endif
460 
461  if (any(y .ne. y)) then
462  call mpistop("NaNs DETECTED IN GCA_INTEGRATE BEFORE ODEINT CALL! ABORTING...")
463  end if
464 
465  ! RK4 integration with adaptive stepwidth
466  call odeint(y,nvar,tloc,tloc+dt_p,eps,h1,hmin,nok,nbad,derivs_gca_rk,rkqs,ierror)
467 
468  if (ierror /= 0) then
469  print *, "odeint returned error code", ierror
470  print *, "1 means hmin too small, 2 means MAXSTP exceeded"
471  print *, "Having a problem with particle", iipart
472  end if
473 
474  if (any(y .ne. y)) then
475  call mpistop("NaNs DETECTED IN GCA_INTEGRATE AFTER ODEINT CALL! ABORTING...")
476  end if
477 
478  ! original RK integration without interpolation in ghost cells
479  ! call odeint(y,nvar,tloc,tloc+dt_p,eps,h1,hmin,nok,nbad,derivs_gca,rkqs)
480 
481  ! final solution vector after rk integration
482  particle(ipart)%self%x(1:ndir) = y(1:ndir)
483  particle(ipart)%self%u(1) = y(ndir+1)
484  particle(ipart)%self%u(2) = y(ndir+2)
485  !particle(ipart)%self%u(3) = y(ndir+3)
486  end select
487 
488  ! now calculate other quantities, mean Lorentz factor, drifts, perpendicular velocity:
489  call get_bfield(igrid_working, y(1:ndir), tloc+dt_p, b)
490  call get_efield(igrid_working, y(1:ndir), tloc+dt_p, e)
491 ! call get_vec(bp, igrid_working,y(1:ndir),tloc+dt_p,b)
492 ! call get_vec(vp, igrid_working,y(1:ndir),tloc+dt_p,vfluid)
493 ! call get_vec(jp, igrid_working,y(1:ndir),tloc+dt_p,current)
494 ! e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
495 ! e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
496 ! e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)
497 
498  absb = sqrt(sum(b(:)**2))
499  bhat(1:ndir) = b(1:ndir) / absb
500 
501  epar = sum(e(:)*bhat(:))
502 
503  call cross(e,bhat,ve)
504 
505  ve(1:ndir) = ve(1:ndir) / absb
506  veabs = sqrt(sum(ve(:)**2))
507  if (relativistic) then
508  kappa = 1.d0/sqrt(1.0d0 - sum(ve(:)**2)/c_norm**2)
509  else
510  kappa = 1.d0
511  end if
512  mr = y(ndir+2); upar = y(ndir+1); m=particle(ipart)%self%m; q=particle(ipart)%self%q
513  if (relativistic) then
514  gamma = sqrt(1.0d0+upar**2/c_norm**2+2.0d0*mr*absb/m/c_norm**2)*kappa
515  else
516  gamma = 1.d0
517  end if
518 
519  particle(ipart)%self%u(3) = gamma
520 
521  ! Time update
522  particle(ipart)%self%time = particle(ipart)%self%time + dt_p
523 
524  ! Update payload
526  particle(ipart)%self%x,particle(ipart)%self%u,q,m,defpayload,ndefpayload,particle(ipart)%self%time)
527  particle(ipart)%payload(1:ndefpayload) = defpayload
528  if (associated(usr_update_payload)) then
530  particle(ipart)%self%x,particle(ipart)%self%u,q,m,usrpayload,nusrpayload,particle(ipart)%self%time)
531  particle(ipart)%payload(ndefpayload+1:npayload) = usrpayload
532  end if
533 
534  end do
535 
536  end subroutine gca_integrate_particles
537 
538  subroutine derivs_gca_rk(t_s,y,dydt)
540 
541  double precision :: t_s, y(ndir+2)
542  double precision :: dydt(ndir+2)
543 
544  double precision,dimension(ndir):: vE, b, e, x, bhat, bdotgradb, vEdotgradb, gradkappaB, vfluid, current
545  double precision,dimension(ndir):: bdotgradvE, vEdotgradvE, u, utmp1, utmp2, utmp3
546  double precision :: upar, Mr, gamma, absb, q, m, epar, kappa
547  integer :: ic^D
548 
549  ! Here the terms in the guiding centre equations of motion are interpolated for
550  ! the RK integration. The interpolation is also done in the ghost cells such
551  ! that the RK integration does not give an error
552 
553  q = particle(ipart_working)%self%q
554  m = particle(ipart_working)%self%m
555 
556  x(1:ndir) = y(1:ndir)
557  upar = y(ndir+1) ! gamma v||
558  mr = y(ndir+2)
559  !gamma = y(ndir+3)
560 
561  if (any(x .ne. x)) then
562  write(*,*) "ERROR IN DERIVS_GCA_RK: NaNs IN X OR Y!"
563  write(*,*) "x",x
564  write(*,*) "y",y(ndir+1:ndir+2)
565  call mpistop("ABORTING...")
566  end if
567 
568  call get_bfield(igrid_working, x, t_s, b)
569  call get_efield(igrid_working, x, t_s, e)
570 ! call get_vec(bp, igrid_working,x,t_s,b)
571 ! call get_vec(vp, igrid_working,x,t_s,vfluid)
572 ! call get_vec(jp, igrid_working,x,t_s,current)
573 ! e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
574 ! e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
575 ! e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)
576  call get_vec(b_dot_grad_b, igrid_working,x,t_s,bdotgradb)
577  call get_vec(ve_dot_grad_b, igrid_working,x,t_s,vedotgradb)
578  call get_vec(grad_kappa_b, igrid_working,x,t_s,gradkappab)
579  call get_vec(b_dot_grad_ve, igrid_working,x,t_s,bdotgradve)
580  call get_vec(ve_dot_grad_ve, igrid_working,x,t_s,vedotgradve)
581 
582  if (any(b .ne. b) .or. any(e .ne. e) &
583  .or. any(bdotgradb .ne. bdotgradb) .or. any(vedotgradb .ne. vedotgradb) &
584  .or. any(gradkappab .ne. gradkappab) .or. any(bdotgradve .ne. bdotgradve) &
585  .or. any(vedotgradve .ne. vedotgradve)) then
586  write(*,*) "ERROR IN DERIVS_GCA_RK: NaNs IN FIELD QUANTITIES!"
587  write(*,*) "b",b
588  write(*,*) "e",e
589  write(*,*) "bdotgradb",bdotgradb
590  write(*,*) "vEdotgradb",vedotgradb
591  write(*,*) "gradkappaB",gradkappab
592  write(*,*) "bdotgradvE",bdotgradve
593  write(*,*) "vEdotgradvE",vedotgradve
594  call mpistop("ABORTING...")
595  end if
596 
597  absb = sqrt(sum(b(:)**2))
598  if (absb .gt. 0.d0) then
599  bhat(1:ndir) = b(1:ndir) / absb
600  else
601  bhat = 0.d0
602  end if
603  epar = sum(e(:)*bhat(:))
604 
605  call cross(e,bhat,ve)
606  if (absb .gt. 0.d0) then
607  ve(1:ndir) = ve(1:ndir) / absb
608  else
609  ve(1:ndir) = 0.d0
610  end if
611 
612  if (relativistic) then
613  kappa = 1.d0/sqrt(1.0d0 - sum(ve(:)**2)/c_norm**2)
614  gamma = sqrt(1.0d0+upar**2/c_norm**2+2.0d0*mr*absb/m/c_norm**2)*kappa
615  else
616  kappa = 1.d0
617  gamma = 1.d0
618  end if
619 
620  if (absb .gt. 0.d0) then
621  utmp1(1:ndir) = bhat(1:ndir)/(absb/kappa**2)
622  else
623  utmp1 = 0.d0
624  end if
625  utmp2(1:ndir) = mr/(gamma*q)*gradkappab(1:ndir) &
626  + m/q* (upar**2/gamma*bdotgradb(1:ndir) + upar*vedotgradb(1:ndir) &
627  + upar*bdotgradve(1:ndir) + gamma*vedotgradve(1:ndir))
628  if (relativistic) then
629  utmp2(1:ndir) = utmp2(1:ndir) + upar*epar/(gamma)*ve(1:ndir)
630  end if
631 
632  call cross(utmp1,utmp2,utmp3)
633  u(1:ndir) = ve(1:ndir) + utmp3(1:ndir)
634 
635  ! done assembling the terms, now write rhs:
636  dydt(1:ndir) = ( u(1:ndir) + upar/gamma * bhat(1:ndir) )
637  dydt(ndir+1) = q/m*epar - mr/(m*gamma) * sum(bhat(:)*gradkappab(:)) &
638  + sum(ve(:)*(upar*bdotgradb(:)+gamma*vedotgradb(:)))
639  dydt(ndir+2) = 0.0d0 ! magnetic moment is conserved
640 
641  end subroutine derivs_gca_rk
642 
643  subroutine derivs_gca(t_s,y,dydt)
645 
646  double precision :: t_s, y(ndir+2)
647  double precision :: dydt(ndir+2)
648 
649  double precision,dimension(ndir):: vE, b, e, x, bhat, bdotgradb, vEdotgradb, gradkappaB, vfluid, current
650  double precision,dimension(ndir):: bdotgradvE, vEdotgradvE, u, utmp1, utmp2, utmp3
651  double precision :: upar, Mr, gamma, absb, q, m, epar, kappa
652 
653  ! Here the normal interpolation is done for the terms in the GCA equations of motion
654 
655  q = particle(ipart_working)%self%q
656  m = particle(ipart_working)%self%m
657 
658  x(1:ndir) = y(1:ndir)
659  upar = y(ndir+1) ! gamma v||
660  mr = y(ndir+2)
661  !gamma = y(ndir+3)
662 
663  if (any(x .ne. x)) then
664  call mpistop("ERROR IN DERIVS_GCA: NaNs IN X! ABORTING...")
665  end if
666 
667  call get_bfield(igrid_working, x, t_s, b)
668  call get_efield(igrid_working, x, t_s, e)
669 ! call get_vec(bp, igrid_working,x,t_s,b)
670 ! call get_vec(vp, igrid_working,x,t_s,vfluid)
671 ! call get_vec(jp, igrid_working,x,t_s,current)
672 ! e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
673 ! e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
674 ! e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)
675  call get_vec(b_dot_grad_b, igrid_working,x,t_s,bdotgradb)
676  call get_vec(ve_dot_grad_b, igrid_working,x,t_s,vedotgradb)
677  call get_vec(grad_kappa_b, igrid_working,x,t_s,gradkappab)
678  call get_vec(b_dot_grad_ve, igrid_working,x,t_s,bdotgradve)
679  call get_vec(ve_dot_grad_ve, igrid_working,x,t_s,vedotgradve)
680 
681  absb = sqrt(sum(b(:)**2))
682  if (absb .gt. 0.d0) then
683  bhat(1:ndir) = b(1:ndir) / absb
684  else
685  bhat = 0.d0
686  end if
687 
688  epar = sum(e(:)*bhat(:))
689  call cross(e,bhat,ve)
690  if (absb .gt. 0.d0) ve(1:ndir) = ve(1:ndir) / absb
691 
692  if (relativistic) then
693  kappa = 1.d0/sqrt(1.0d0 - sum(ve(:)**2)/c_norm**2)
694  gamma = sqrt(1.0d0+upar**2/c_norm**2+2.0d0*mr*absb/m/c_norm**2)*kappa
695  else
696  kappa = 1.d0
697  gamma = 1.d0
698  end if
699  if (absb .gt. 0.d0) then
700  utmp1(1:ndir) = bhat(1:ndir)/(absb/kappa**2)
701  else
702  utmp1 = 0.d0
703  end if
704  utmp2(1:ndir) = mr/(gamma*q)*gradkappab(1:ndir) &
705  + m/q* (upar**2/gamma*bdotgradb(1:ndir) + upar*vedotgradb(1:ndir) &
706  + upar*bdotgradve(1:ndir) + gamma*vedotgradve(1:ndir))
707  if (relativistic) then
708  utmp2(1:ndir) = utmp2(1:ndir) + upar*epar/(gamma)*ve(1:ndir)
709  end if
710 
711  call cross(utmp1,utmp2,utmp3)
712  u(1:ndir) = ve(1:ndir) + utmp3(1:ndir)
713 
714  ! done assembling the terms, now write rhs:
715  dydt(1:ndir) = ( u(1:ndir) + upar/gamma * bhat(1:ndir) )
716  dydt(ndir+1) = q/m*epar - mr/(m*gamma) * sum(bhat(:)*gradkappab(:)) &
717  + sum(ve(:)*(upar*bdotgradb(:)+gamma*vedotgradb(:)))
718  dydt(ndir+2) = 0.0d0 ! magnetic moment is conserved
719 
720  end subroutine derivs_gca
721 
722  !> Update payload subroutine
723  subroutine gca_update_payload(igrid,xpart,upart,qpart,mpart,mypayload,mynpayload,particle_time)
725  integer, intent(in) :: igrid,mynpayload
726  double precision, intent(in) :: xpart(1:ndir),upart(1:ndir),qpart,mpart,particle_time
727  double precision, intent(out) :: mypayload(mynpayload)
728  double precision, dimension(1:ndir) :: vE, e, b, bhat, vfluid, current
729  double precision, dimension(1:ndir) :: drift1, drift2
730  double precision, dimension(1:ndir) :: drift3, drift4, drift5, drift6, drift7
731  double precision, dimension(1:ndir) :: bdotgradb, vEdotgradb, gradkappaB
732  double precision, dimension(1:ndir) :: bdotgradvE, vEdotgradvE
733  double precision, dimension(1:ndir) :: gradBdrift, reldrift, bdotgradbdrift
734  double precision, dimension(1:ndir) :: vEdotgradbdrift, bdotgradvEdrift
735  double precision, dimension(1:ndir) :: vEdotgradvEdrift
736  double precision :: kappa, upar, absb, gamma, vpar, vEabs
737  double precision :: gradBdrift_abs, reldrift_abs, epar
738  double precision :: bdotgradbdrift_abs, vEdotgradbdrift_abs
739  double precision :: bdotgradvEdrift_abs, vEdotgradvEdrift_abs
740  double precision :: momentumpar1, momentumpar2, momentumpar3, momentumpar4
741 
742  call get_bfield(igrid, xpart(1:ndir), particle_time, b)
743  call get_efield(igrid, xpart(1:ndir), particle_time, e)
744 ! call get_vec(bp, igrid,xpart(1:ndir),particle_time,b)
745 ! call get_vec(vp, igrid,xpart(1:ndir),particle_time,vfluid)
746 ! call get_vec(jp, igrid,xpart(1:ndir),particle_time,current)
747 ! e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
748 ! e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
749 ! e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)
750 
751  absb = sqrt(sum(b(:)**2))
752  bhat(1:ndir) = b(1:ndir) / absb
753  epar = sum(e(:)*bhat(:))
754  call cross(e,bhat,ve)
755  ve(1:ndir) = ve(1:ndir) / absb
756  veabs = sqrt(sum(ve(:)**2))
757  if (relativistic) then
758  kappa = 1.d0/sqrt(1.0d0 - sum(ve(:)**2)/c_norm**2)
759  else
760  kappa = 1.d0
761  end if
762  vpar = upart(1)/upart(3)
763  upar = upart(1)
764 
765  call get_vec(b_dot_grad_b, igrid,xpart(1:ndir),particle_time,bdotgradb)
766  call get_vec(ve_dot_grad_b, igrid,xpart(1:ndir),particle_time,vedotgradb)
767  call get_vec(grad_kappa_b, igrid,xpart(1:ndir),particle_time,gradkappab)
768  call get_vec(b_dot_grad_ve, igrid,xpart(1:ndir),particle_time,bdotgradve)
769  call get_vec(ve_dot_grad_ve, igrid,xpart(1:ndir),particle_time,vedotgradve)
770 
771  drift1(1:ndir) = bhat(1:ndir)/(absb/kappa**2)
772  drift2(1:ndir) = upart(2)/(upart(3)*qpart)*gradkappab(1:ndir)
773 
774  call cross(drift1,drift2,gradbdrift)
775  gradbdrift_abs = sqrt(sum(gradbdrift(:)**2))
776 
777  drift3(1:ndir) = upar*epar/upart(3)*ve(1:ndir)
778  call cross(drift1,drift3,reldrift)
779  reldrift_abs = sqrt(sum(reldrift(:)**2))
780 
781  drift4(1:ndir) = mpart/qpart* ( upar**2/upart(3)*bdotgradb(1:ndir))
782  call cross(drift1,drift4,bdotgradbdrift)
783  bdotgradbdrift_abs = sqrt(sum(bdotgradbdrift(:)**2))
784 
785  drift5(1:ndir) = mpart/qpart* ( upar*vedotgradb(1:ndir))
786  call cross(drift1,drift5,vedotgradbdrift)
787  vedotgradbdrift_abs = sqrt(sum(vedotgradbdrift(:)**2))
788 
789  drift6(1:ndir) = mpart/qpart* ( upar*bdotgradve(1:ndir))
790  call cross(drift1,drift6,bdotgradvedrift)
791  bdotgradvedrift_abs = sqrt(sum(bdotgradvedrift(:)**2))
792 
793  drift7(1:ndir) = mpart/qpart* (upart(3)*vedotgradve(1:ndir))
794  call cross(drift1,drift7,vedotgradvedrift)
795  vedotgradvedrift_abs = sqrt(sum(vedotgradvedrift(:)**2))
796 
797  momentumpar1 = qpart/mpart*epar
798  momentumpar2 = -(upart(2)/mpart/upart(3))*sum(bhat(:)*gradkappab(:))
799  momentumpar3 = upar*sum(ve(:)*bdotgradb(:))
800  momentumpar4 = upart(3)*sum(ve(:)*vedotgradb(:))
801 
802  ! Payload update
803  if (mynpayload > 0) then
804  ! current gyroradius
805  mypayload(1) = sqrt(2.0d0*mpart*upart(2)*absb)/abs(qpart*absb)
806  end if
807  if (mynpayload > 1) then
808  ! pitch angle
809  mypayload(2) = atan2(sqrt((2.0d0*upart(2)*absb)/(mpart*upart(3)**2)),vpar)
810  end if
811  if (mynpayload > 2) then
812  ! particle v_perp
813  mypayload(3) = sqrt((2.0d0*upart(2)*absb)/(mpart*upart(3)**2))
814  end if
815  if (mynpayload > 3) then
816  ! particle parallel momentum term 1
817  mypayload(4) = momentumpar1
818  end if
819  if (mynpayload > 4) then
820  ! particle parallel momentum term 2
821  mypayload(5) = momentumpar2
822  end if
823  if (mynpayload > 5) then
824  ! particle parallel momentum term 3
825  mypayload(6) = momentumpar3
826  end if
827  if (mynpayload > 6) then
828  ! particle parallel momentum term 4
829  mypayload(7) = momentumpar4
830  end if
831  if (mynpayload > 7) then
832  ! particle ExB drift
833  mypayload(8) = veabs
834  end if
835  if (mynpayload > 8) then
836  ! relativistic drift
837  mypayload(9) = gradbdrift_abs
838  end if
839  if (mynpayload > 9) then
840  ! gradB drift
841  mypayload(10) = reldrift_abs
842  end if
843  if (mynpayload > 10) then
844  ! bdotgradb drift
845  mypayload(11) = bdotgradbdrift_abs
846  end if
847  if (mynpayload > 11) then
848  ! vEdotgradb drift
849  mypayload(12) = vedotgradbdrift_abs
850  end if
851  if (mynpayload > 12) then
852  ! bdotgradvE drift
853  mypayload(13) = bdotgradvedrift_abs
854  end if
855  if (mynpayload > 13) then
856  ! vEdotgradvE drift
857  mypayload(14) = vedotgradvedrift_abs
858  end if
859 
860  end subroutine gca_update_payload
861 
862  function gca_get_particle_dt(partp, end_time) result(dt_p)
863  use mod_odeint
865  type(particle_ptr), intent(in) :: partp
866  double precision, intent(in) :: end_time
867  double precision :: dt_p
868 
869  double precision :: tout, dt_particles_mype, dt_cfl0, dt_cfl1, dt_a
870  double precision :: dxmin, vp, a, gammap
871  double precision :: v(ndir), y(ndir+2),ytmp(ndir+2), dydt(ndir+2), v0(ndir), v1(ndir), dydt1(ndir+2)
872  double precision :: ap0, ap1, dt_cfl_ap0, dt_cfl_ap1, dt_cfl_ap
873  double precision :: dt_euler, dt_tmp
874  ! make these particle cfl conditions more restrictive if you are interpolating out of the grid
875  double precision :: cfl, uparcfl
876  double precision, parameter :: uparmin=1.0d-6*const_c
877  integer :: ipart, iipart, nout, ic^d, igrid_particle, ipe_particle, ipe
878  logical :: bc_applied
879 
880  if (const_dt_particles > 0) then
881  dt_p = const_dt_particles
882  return
883  end if
884 
885  cfl = particles_cfl
886  uparcfl = particles_cfl
887 
888  igrid_working = partp%igrid
889  ipart_working = partp%self%index
890  dt_tmp = (end_time - partp%self%time)
891  if(dt_tmp .le. 0.0d0) dt_tmp = smalldouble
892  ! make sure we step only one cell at a time, first check CFL at current location
893  ! then we make an Euler step to the new location and check the new CFL
894  ! we simply take the minimum of the two timesteps.
895  ! added safety factor cfl:
896  dxmin = min({rnode(rpdx^d_,igrid_working)},bigdouble)*cfl
897  ! initial solution vector:
898  y(1:ndir) = partp%self%x(1:ndir) ! position of guiding center
899  y(ndir+1) = partp%self%u(1) ! parallel momentum component (gamma v||)
900  y(ndir+2) = partp%self%u(2) ! conserved magnetic moment
901  ytmp=y
902  !y(ndir+3) = partp%self%u(3) ! Lorentz factor of guiding centre
903 
904  call derivs_gca(partp%self%time,y,dydt)
905  v0(1:ndir) = dydt(1:ndir)
906  ap0 = dydt(ndir+1)
907 
908  ! guiding center velocity:
909  v(1:ndir) = abs(dydt(1:ndir))
910  vp = sqrt(sum(v(:)**2))
911 
912  dt_cfl0 = dxmin / max(vp, smalldouble)
913  dt_cfl_ap0 = uparcfl * abs(max(abs(y(ndir+1)),uparmin) / max(abs(ap0), smalldouble))
914  !dt_cfl_ap0 = min(dt_cfl_ap0, uparcfl * sqrt(abs(unit_length*dxmin/(ap0+smalldouble))) )
915 
916  ! make an Euler step with the proposed timestep:
917  ! new solution vector:
918  dt_euler = min(dt_tmp,dt_cfl0,dt_cfl_ap0)
919  y(1:ndir+2) = y(1:ndir+2) + dt_euler * dydt(1:ndir+2)
920 
921  partp%self%x(1:ndir) = y(1:ndir) ! position of guiding center
922  partp%self%u(1) = y(ndir+1) ! parallel momentum component (gamma v||)
923  partp%self%u(2) = y(ndir+2) ! conserved magnetic moment
924 
925  ! first check if the particle is outside the physical domain or in the ghost cells
927  y(1:ndir+2) = ytmp(1:ndir+2)
928  end if
929 
930  call derivs_gca_rk(partp%self%time+dt_euler,y,dydt)
931  !call derivs_gca(partp%self%time+dt_euler,y,dydt)
932 
933  v1(1:ndir) = dydt(1:ndir)
934  ap1 = dydt(ndir+1)
935 
936  ! guiding center velocity:
937  v(1:ndir) = abs(dydt(1:ndir))
938  vp = sqrt(sum(v(:)**2))
939 
940  dt_cfl1 = dxmin / max(vp, smalldouble)
941  dt_cfl_ap1 = uparcfl * abs(max(abs(y(ndir+1)),uparmin) / max(abs(ap1), smalldouble))
942  !dt_cfl_ap1 = min(dt_cfl_ap1, uparcfl * sqrt(abs(unit_length*dxmin/(ap1+smalldouble))) )
943 
944  dt_tmp = min(dt_euler, dt_cfl1, dt_cfl_ap1)
945 
946  partp%self%x(1:ndir) = ytmp(1:ndir) ! position of guiding center
947  partp%self%u(1) = ytmp(ndir+1) ! parallel momentum component (gamma v||)
948  partp%self%u(2) = ytmp(ndir+2) ! conserved magnetic moment
949  !dt_tmp = min(dt_cfl1, dt_cfl_ap1)
950 
951  ! time step due to parallel acceleration:
952  ! The standard thing, dt=sqrt(dx/a) where we compute a from d(gamma v||)/dt and d(gamma)/dt
953  ! dt_ap = sqrt(abs(dxmin*unit_length*y(ndir+3)/( dydt(ndir+1) - y(ndir+1)/y(ndir+3)*dydt(ndir+3) ) ) )
954  ! vp = sqrt(sum(v(1:ndir)**))
955  ! gammap = sqrt(1.0d0/(1.0d0-(vp/const_c)**2))
956  ! ap = const_c**2/vp*gammap**(-3)*dydt(ndir+3)
957  ! dt_ap = sqrt(dxmin*unit_length/ap)
958 
959  !dt_a = bigdouble
960  !if (dt_euler .gt. smalldouble) then
961  ! a = sqrt(sum((v1(1:ndir)-v0(1:ndir))**2))/dt_euler
962  ! dt_a = min(sqrt(dxmin/a),bigdouble)
963  !end if
964 
965  !dt_p = min(dt_tmp , dt_a)
966  dt_p = dt_tmp
967 
968  ! Make sure we do not advance beyond end_time
969  call limit_dt_endtime(end_time - partp%self%time, dt_p)
970 
971  end function gca_get_particle_dt
972 
973 end module mod_particle_gca
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
Definition: mod_geometry.t:321
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 icomm
The MPI communicator.
integer mype
The rank of the current MPI task.
integer, 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.
Definition: mod_odeint.t:2
subroutine odeint(ystart, nvar, x1, x2, eps, h1, hmin, nok, nbad, derivs, rkqs, ierror)
Definition: mod_odeint.t:13
subroutine rkqs(y, dydx, n, x, htry, eps, yscal, hdid, hnext, derivs)
Definition: mod_odeint.t:116
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.
logical function particle_in_igrid(ipart, igrid)
Quick check if particle is still in igrid.
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
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 derivs_gca(t_s, y, dydt)
subroutine derivs_gca_rk(t_s, y, dydt)
subroutine, public gca_init()
double precision function gca_get_particle_dt(partp, end_time)
subroutine, public gca_create_particles()
subroutine gca_fill_gridvars
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.
subroutine gca_integrate_particles(end_time)
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...
subroutine gca_update_payload(igrid, xpart, upart, qpart, mpart, mypayload, mynpayload, particle_time)
Update payload subroutine.
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