MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
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
31contains
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
91 particles_fill_gridvars => gca_fill_gridvars
92
93 if (associated(particles_define_additional_gridvars)) then
95 end if
96
97 select case(integrator_type_particles)
98 case('RK4','Rk4','rk4')
100 case('ARK4','ARk4','Ark4','ark4')
101 integrator = ark4
102 case default
103 integrator = ark4
104 end select
105
106 particles_integrate => gca_integrate_particles
107
108 end subroutine gca_init
109
111 ! initialise the particles
114
115 double precision :: b(ndir), u(ndir), magmom
116 double precision :: bnorm, lfac, vnorm, vperp, vpar
117 double precision :: x(3, num_particles)
118 double precision :: v(3, num_particles)
119 double precision :: q(num_particles)
120 double precision :: m(num_particles)
121 double precision :: rrd(num_particles,ndir)
122 double precision :: defpayload(ndefpayload)
123 double precision :: usrpayload(nusrpayload)
124 integer :: igrid, ipe_particle
125 integer :: n, idir, nparticles_local
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 ! write the first csv file of particles
219
220 end subroutine gca_create_particles
221
222 subroutine gca_fill_gridvars
225 use mod_geometry
226
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
233 integer :: idirmin
234
236
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)
240
241 ! grad(kappa B)
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))
249 do idir=1,ndir
250 where (absb(ixg^t) .gt. 0.d0)
251 ve(ixg^t,idir) = ve(ixg^t,idir) / absb(ixg^t)**2
252 elsewhere
253 ve(ixg^t,idir) = 0.d0
254 end where
255 end do
256 if (any(sum(ve(ixg^t,:)**2,dim=ndim+1) .ge. c_norm**2) .and. relativistic) then
257 call mpistop("GCA FILL GRIDVARS: vE>c! ABORTING...")
258 end if
259 if (any(ve .ne. ve)) then
260 call mpistop("GCA FILL GRIDVARS: NaNs IN vE! ABORTING...")
261 end if
262
263 if (relativistic) then
264 kappa(ixg^t) = 1.d0/sqrt(1.0d0 - sum(ve(ixg^t,:)**2,dim=ndim+1)/c_norm**2)
265 else
266 kappa(ixg^t) = 1.d0
267 end if
268 kappa_b(ixg^t) = absb(ixg^t) / kappa(ixg^t)
269
270 if (any(kappa_b .ne. kappa_b)) then
271 call mpistop("GCA FILL GRIDVARS: NaNs IN kappa_B! ABORTING...")
272 end if
273
274 tmp=0.d0
275 do idim=1,ndim
276 call gradient(kappa_b,ixg^ll,ixg^ll^lsub1,idim,tmp)
277 gridvars(igrid)%w(ixg^t,grad_kappa_b(idim)) = tmp(ixg^t)
278 end do
279
280 ! bhat
281 do idir=1,ndir
282 where (absb(ixg^t) .gt. 0.d0)
283 bhat(ixg^t,idir) = gridvars(igrid)%w(ixg^t,bp(idir)) / absb(ixg^t)
284 elsewhere
285 bhat(ixg^t,idir) = 0.d0
286 end where
287 end do
288 if (any(bhat .ne. bhat)) then
289 call mpistop("GCA FILL GRIDVARS: NaNs IN bhat! ABORTING...")
290 end if
291
292 do idir=1,ndir
293 ! (b dot grad) b and the other directional derivatives
294 do idim=1,ndim
295 call gradient(bhat(ixg^t,idir),ixg^ll,ixg^ll^lsub1,idim,tmp)
296 gridvars(igrid)%w(ixg^t,b_dot_grad_b(idir)) = gridvars(igrid)%w(ixg^t,b_dot_grad_b(idir)) &
297 + bhat(ixg^t,idim) * tmp(ixg^t)
298 gridvars(igrid)%w(ixg^t,ve_dot_grad_b(idir)) = gridvars(igrid)%w(ixg^t,ve_dot_grad_b(idir)) &
299 + ve(ixg^t,idim) * tmp(ixg^t)
300 call gradient(ve(ixg^t,idir),ixg^ll,ixg^ll^lsub1,idim,tmp)
301 gridvars(igrid)%w(ixg^t,b_dot_grad_ve(idir)) = gridvars(igrid)%w(ixg^t,b_dot_grad_ve(idir)) &
302 + bhat(ixg^t,idim) * tmp(ixg^t)
303 gridvars(igrid)%w(ixg^t,ve_dot_grad_ve(idir)) = gridvars(igrid)%w(ixg^t,ve_dot_grad_ve(idir)) &
304 + ve(ixg^t,idim) * tmp(ixg^t)
305 end do
306 end do
307
308 end do
309
310 end subroutine gca_fill_gridvars
311
312 subroutine gca_integrate_particles(end_time)
313 use mod_odeint
316 double precision, intent(in) :: end_time
317
318 double precision :: lfac, abss
319 double precision :: defpayload(ndefpayload)
320 double precision :: usrpayload(nusrpayload)
321 double precision :: dt_p, tloc, y(ndir+2),dydt(ndir+2),ytmp(ndir+2), euler_cfl, int_factor
322 double precision :: tk, k1(ndir+2),k2(ndir+2),k3(ndir+2),k4(ndir+2)
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
336 ! Precision of time-integration:
337 double precision,parameter :: eps=1.0d-6
338 ! for odeint:
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
344
345 nvar=ndir+2
346
347 do iipart=1,nparticles_active_on_mype
348 ipart = particles_active_on_mype(iipart)
349 int_choice = .false.
350 dt_p = gca_get_particle_dt(particle(ipart), end_time)
351 particle(ipart)%self%dt = dt_p
352
353 igrid_working = particle(ipart)%igrid
354 ipart_working = particle(ipart)%self%index
355 tloc = particle(ipart)%self%time
356 x(1:ndir) = particle(ipart)%self%x(1:ndir)
357
358 select case (integrator)
359 case (rk4)
360 ! Runge-Kutta order 4
361 tk = tloc
362 y(1:ndir) = x(1:ndir) ! position of guiding center
363 y(ndir+1) = particle(ipart)%self%u(1) ! parallel momentum component (gamma v||)
364 y(ndir+2) = particle(ipart)%self%u(2) ! conserved magnetic moment Mr
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
368 y(ndir+1) = particle(ipart)%self%u(1) + dt_p*k1(ndir+1)/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
372 y(ndir+1) = particle(ipart)%self%u(1) + dt_p*k2(ndir+1)/2.d0
373 call derivs_gca_rk(tk,y,k3)
374 tk = tloc + dt_p
375 y(1:ndir) = x + dt_p*k3(1:ndir)
376 y(ndir+1) = particle(ipart)%self%u(1) + dt_p*k3(ndir+1)
377 call derivs_gca_rk(tk,y,k4)
378 y(1:ndir) = x(1:ndir) ! position of guiding center
379 y(ndir+1) = particle(ipart)%self%u(1) ! parallel momentum component (gamma v||)
380 y = y + dt_p/6.d0*(k1 + 2.d0*k2 + 2.d0*k3 + k4)
381 particle(ipart)%self%x(1:ndir) = y(1:ndir)
382 particle(ipart)%self%u(1) = y(ndir+1)
383
384 case (ark4)
385 ! Adaptive stepwidth RK4:
386 ! initial solution vector:
387 y(1:ndir) = x(1:ndir) ! position of guiding center
388 y(ndir+1) = particle(ipart)%self%u(1) ! parallel momentum component (gamma v||)
389 y(ndir+2) = particle(ipart)%self%u(2) ! conserved magnetic moment Mr
390 ! y(ndir+3) = particle(ipart)%self%u(3) ! Lorentz factor of particle
391
392 ! we temporarily save the solution vector, to replace the one from the euler
393 ! timestep after euler integration
394 ytmp=y
395
396 call derivs_gca(particle(ipart)%self%time,y,dydt)
397
398 ! make an Euler step with the proposed timestep:
399 ! factor to ensure we capture all particles near the internal ghost cells.
400 ! Can be adjusted during a run, after an interpolation error.
401 euler_cfl=2.5d0
402
403 ! new solution vector:
404 y(1:ndir+2) = y(1:ndir+2) + euler_cfl * dt_p * dydt(1:ndir+2)
405 particle(ipart)%self%x(1:ndir) = y(1:ndir) ! position of guiding center
406 particle(ipart)%self%u(1) = y(ndir+1) ! parallel momentum component(gamma v||)
407 particle(ipart)%self%u(2) = y(ndir+2) ! conserved magnetic moment
408
409 ! check if the particle is in the internal ghost cells
410 int_factor =1.0d0
411
413 ! if particle is not in the grid an euler timestep is taken instead of a RK4
414 ! timestep. Then based on that we do an interpolation and check how much further
415 ! the timestep for the RK4 has to be restricted.
416 ! factor to make integration more accurate for particles near the internal
417 ! ghost cells. This factor can be changed during integration after an
418 ! interpolation error. But one should be careful with timesteps for i/o
419
420 ! flat interpolation:
421 {ic^d = int((y(^d)-rnode(rpxmin^d_,igrid_working))/rnode(rpdx^d_,igrid_working)) + 1 + nghostcells\}
422
423 ! linear interpolation:
424 {
425 if (ps(igrid_working)%x({ic^dd},^d) .lt. y(^d)) then
426 ic1^d = ic^d
427 else
428 ic1^d = ic^d -1
429 end if
430 ic2^d = ic1^d + 1
431 \}
432
433 int_factor =0.5d0
434
435 {^d&
436 if (ic1^d .le. ixglo^d-2 .or. ic2^d .ge. ixghi^d+2) then
437 int_factor = 0.05d0
438 end if
439 \}
440
441 {^d&
442 if (ic1^d .eq. ixglo^d-1 .or. ic2^d .eq. ixghi^d+1) then
443 int_factor = 0.1d0
444 end if
445 \}
446
447 dt_p=int_factor*dt_p
448 end if
449
450 ! replace the solution vector with the original as it was before the Euler timestep
451 y(1:ndir+2) = ytmp(1:ndir+2)
452
453 particle(ipart)%self%x(1:ndir) = ytmp(1:ndir) ! position of guiding center
454 particle(ipart)%self%u(1) = ytmp(ndir+1) ! parallel momentum component (gamma v||)
455 particle(ipart)%self%u(2) = ytmp(ndir+2) ! conserved magnetic moment
456
457 ! specify a minimum step hmin. If the timestep reaches this minimum, multiply by
458 ! a factor 100 to make sure the RK integration doesn't crash
459 h1 = dt_p/2.0d0; hmin=1.0d-9; h_old=dt_p/2.0d0
460
461 if(h1 .lt. hmin)then
462 h1=hmin
463 dt_p=2.0d0*h1
464 endif
465
466 if (any(y .ne. y)) then
467 call mpistop("NaNs DETECTED IN GCA_INTEGRATE BEFORE ODEINT CALL! ABORTING...")
468 end if
469
470 ! RK4 integration with adaptive stepwidth
471 call odeint(y,nvar,tloc,tloc+dt_p,eps,h1,hmin,nok,nbad,derivs_gca_rk,rkqs,ierror)
472
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
477 end if
478
479 if (any(y .ne. y)) then
480 call mpistop("NaNs DETECTED IN GCA_INTEGRATE AFTER ODEINT CALL! ABORTING...")
481 end if
482
483 ! original RK integration without interpolation in ghost cells
484 ! call odeint(y,nvar,tloc,tloc+dt_p,eps,h1,hmin,nok,nbad,derivs_gca,rkqs)
485
486 ! final solution vector after rk integration
487 particle(ipart)%self%x(1:ndir) = y(1:ndir)
488 particle(ipart)%self%u(1) = y(ndir+1)
489 particle(ipart)%self%u(2) = y(ndir+2)
490 !particle(ipart)%self%u(3) = y(ndir+3)
491 end select
492
493 ! now calculate other quantities, mean Lorentz factor, drifts, perpendicular velocity:
494 call get_bfield(igrid_working, y(1:ndir), tloc+dt_p, b)
495 call get_efield(igrid_working, y(1:ndir), tloc+dt_p, e)
496! call get_vec(bp, igrid_working,y(1:ndir),tloc+dt_p,b)
497! call get_vec(vp, igrid_working,y(1:ndir),tloc+dt_p,vfluid)
498! call get_vec(jp, igrid_working,y(1:ndir),tloc+dt_p,current)
499! e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
500! e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
501! e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)
502
503 absb = sqrt(sum(b(:)**2))
504 bhat(1:ndir) = b(1:ndir) / absb
505
506 epar = sum(e(:)*bhat(:))
507
508 call cross(e,bhat,ve)
509
510 ve(1:ndir) = ve(1:ndir) / absb
511 veabs = sqrt(sum(ve(:)**2))
512 if (relativistic) then
513 kappa = 1.d0/sqrt(1.0d0 - sum(ve(:)**2)/c_norm**2)
514 else
515 kappa = 1.d0
516 end if
517 mr = y(ndir+2); upar = y(ndir+1); m=particle(ipart)%self%m; q=particle(ipart)%self%q
518 if (relativistic) then
519 gamma = sqrt(1.0d0+upar**2/c_norm**2+2.0d0*mr*absb/m/c_norm**2)*kappa
520 else
521 gamma = 1.d0
522 end if
523
524 particle(ipart)%self%u(3) = gamma
525
526 ! Time update
527 particle(ipart)%self%time = particle(ipart)%self%time + dt_p
528
529 ! Update payload
530 call gca_update_payload(igrid_working,&
531 particle(ipart)%self%x,particle(ipart)%self%u,q,m,defpayload,ndefpayload,particle(ipart)%self%time)
532 particle(ipart)%payload(1:ndefpayload) = defpayload
533 if (associated(usr_update_payload)) then
535 particle(ipart)%self%x,particle(ipart)%self%u,q,m,usrpayload,nusrpayload,particle(ipart)%self%time)
536 particle(ipart)%payload(ndefpayload+1:npayload) = usrpayload
537 end if
538
539 end do
540
541 end subroutine gca_integrate_particles
542
543 subroutine derivs_gca_rk(t_s,y,dydt)
545
546 double precision :: t_s, y(ndir+2)
547 double precision :: dydt(ndir+2)
548
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
552 integer :: ic^d
553
554 ! Here the terms in the guiding centre equations of motion are interpolated for
555 ! the RK integration. The interpolation is also done in the ghost cells such
556 ! that the RK integration does not give an error
557
558 q = particle(ipart_working)%self%q
559 m = particle(ipart_working)%self%m
560
561 x(1:ndir) = y(1:ndir)
562 upar = y(ndir+1) ! gamma v||
563 mr = y(ndir+2)
564 !gamma = y(ndir+3)
565
566 if (any(x .ne. x)) then
567 write(*,*) "ERROR IN DERIVS_GCA_RK: NaNs IN X OR Y!"
568 write(*,*) "x",x
569 write(*,*) "y",y(ndir+1:ndir+2)
570 call mpistop("ABORTING...")
571 end if
572
573 call get_bfield(igrid_working, x, t_s, b)
574 call get_efield(igrid_working, x, t_s, e)
575! call get_vec(bp, igrid_working,x,t_s,b)
576! call get_vec(vp, igrid_working,x,t_s,vfluid)
577! call get_vec(jp, igrid_working,x,t_s,current)
578! e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
579! e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
580! e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)
581 call get_vec(b_dot_grad_b, igrid_working,x,t_s,bdotgradb)
582 call get_vec(ve_dot_grad_b, igrid_working,x,t_s,vedotgradb)
583 call get_vec(grad_kappa_b, igrid_working,x,t_s,gradkappab)
584 call get_vec(b_dot_grad_ve, igrid_working,x,t_s,bdotgradve)
585 call get_vec(ve_dot_grad_ve, igrid_working,x,t_s,vedotgradve)
586
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!"
592 write(*,*) "b",b
593 write(*,*) "e",e
594 write(*,*) "bdotgradb",bdotgradb
595 write(*,*) "vEdotgradb",vedotgradb
596 write(*,*) "gradkappaB",gradkappab
597 write(*,*) "bdotgradvE",bdotgradve
598 write(*,*) "vEdotgradvE",vedotgradve
599 call mpistop("ABORTING...")
600 end if
601
602 absb = sqrt(sum(b(:)**2))
603 if (absb .gt. 0.d0) then
604 bhat(1:ndir) = b(1:ndir) / absb
605 else
606 bhat = 0.d0
607 end if
608 epar = sum(e(:)*bhat(:))
609
610 call cross(e,bhat,ve)
611 if (absb .gt. 0.d0) then
612 ve(1:ndir) = ve(1:ndir) / absb
613 else
614 ve(1:ndir) = 0.d0
615 end if
616
617 if (relativistic) 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
620 else
621 kappa = 1.d0
622 gamma = 1.d0
623 end if
624
625 if (absb .gt. 0.d0) then
626 utmp1(1:ndir) = bhat(1:ndir)/(absb/kappa**2)
627 else
628 utmp1 = 0.d0
629 end if
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))
633 if (relativistic) then
634 utmp2(1:ndir) = utmp2(1:ndir) + upar*epar/(gamma)*ve(1:ndir)
635 end if
636
637 call cross(utmp1,utmp2,utmp3)
638 u(1:ndir) = ve(1:ndir) + utmp3(1:ndir)
639
640 ! done assembling the terms, now write rhs:
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(:)))
644 dydt(ndir+2) = 0.0d0 ! magnetic moment is conserved
645
646 end subroutine derivs_gca_rk
647
648 subroutine derivs_gca(t_s,y,dydt)
650
651 double precision :: t_s, y(ndir+2)
652 double precision :: dydt(ndir+2)
653
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
657
658 ! Here the normal interpolation is done for the terms in the GCA equations of motion
659
660 q = particle(ipart_working)%self%q
661 m = particle(ipart_working)%self%m
662
663 x(1:ndir) = y(1:ndir)
664 upar = y(ndir+1) ! gamma v||
665 mr = y(ndir+2)
666 !gamma = y(ndir+3)
667
668 if (any(x .ne. x)) then
669 call mpistop("ERROR IN DERIVS_GCA: NaNs IN X! ABORTING...")
670 end if
671
672 call get_bfield(igrid_working, x, t_s, b)
673 call get_efield(igrid_working, x, t_s, e)
674! call get_vec(bp, igrid_working,x,t_s,b)
675! call get_vec(vp, igrid_working,x,t_s,vfluid)
676! call get_vec(jp, igrid_working,x,t_s,current)
677! e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
678! e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
679! e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)
680 call get_vec(b_dot_grad_b, igrid_working,x,t_s,bdotgradb)
681 call get_vec(ve_dot_grad_b, igrid_working,x,t_s,vedotgradb)
682 call get_vec(grad_kappa_b, igrid_working,x,t_s,gradkappab)
683 call get_vec(b_dot_grad_ve, igrid_working,x,t_s,bdotgradve)
684 call get_vec(ve_dot_grad_ve, igrid_working,x,t_s,vedotgradve)
685
686 absb = sqrt(sum(b(:)**2))
687 if (absb .gt. 0.d0) then
688 bhat(1:ndir) = b(1:ndir) / absb
689 else
690 bhat = 0.d0
691 end if
692
693 epar = sum(e(:)*bhat(:))
694 call cross(e,bhat,ve)
695 if (absb .gt. 0.d0) ve(1:ndir) = ve(1:ndir) / absb
696
697 if (relativistic) then
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
700 else
701 kappa = 1.d0
702 gamma = 1.d0
703 end if
704 if (absb .gt. 0.d0) then
705 utmp1(1:ndir) = bhat(1:ndir)/(absb/kappa**2)
706 else
707 utmp1 = 0.d0
708 end if
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))
712 if (relativistic) then
713 utmp2(1:ndir) = utmp2(1:ndir) + upar*epar/(gamma)*ve(1:ndir)
714 end if
715
716 call cross(utmp1,utmp2,utmp3)
717 u(1:ndir) = ve(1:ndir) + utmp3(1:ndir)
718
719 ! done assembling the terms, now write rhs:
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(:)))
723 dydt(ndir+2) = 0.0d0 ! magnetic moment is conserved
724
725 end subroutine derivs_gca
726
727 !> Update payload subroutine
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
746
747 call get_bfield(igrid, xpart(1:ndir), particle_time, b)
748 call get_efield(igrid, xpart(1:ndir), particle_time, e)
749! call get_vec(bp, igrid,xpart(1:ndir),particle_time,b)
750! call get_vec(vp, igrid,xpart(1:ndir),particle_time,vfluid)
751! call get_vec(jp, igrid,xpart(1:ndir),particle_time,current)
752! e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
753! e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
754! e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)
755
756 absb = sqrt(sum(b(:)**2))
757 bhat(1:ndir) = b(1:ndir) / absb
758 epar = sum(e(:)*bhat(:))
759 call cross(e,bhat,ve)
760 ve(1:ndir) = ve(1:ndir) / absb
761 veabs = sqrt(sum(ve(:)**2))
762 if (relativistic) then
763 kappa = 1.d0/sqrt(1.0d0 - sum(ve(:)**2)/c_norm**2)
764 else
765 kappa = 1.d0
766 end if
767 vpar = upart(1)/upart(3)
768 upar = upart(1)
769
770 call get_vec(b_dot_grad_b, igrid,xpart(1:ndir),particle_time,bdotgradb)
771 call get_vec(ve_dot_grad_b, igrid,xpart(1:ndir),particle_time,vedotgradb)
772 call get_vec(grad_kappa_b, igrid,xpart(1:ndir),particle_time,gradkappab)
773 call get_vec(b_dot_grad_ve, igrid,xpart(1:ndir),particle_time,bdotgradve)
774 call get_vec(ve_dot_grad_ve, igrid,xpart(1:ndir),particle_time,vedotgradve)
775
776 drift1(1:ndir) = bhat(1:ndir)/(absb/kappa**2)
777 drift2(1:ndir) = upart(2)/(upart(3)*qpart)*gradkappab(1:ndir)
778
779 call cross(drift1,drift2,gradbdrift)
780 gradbdrift_abs = sqrt(sum(gradbdrift(:)**2))
781
782 drift3(1:ndir) = upar*epar/upart(3)*ve(1:ndir)
783 call cross(drift1,drift3,reldrift)
784 reldrift_abs = sqrt(sum(reldrift(:)**2))
785
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))
789
790 drift5(1:ndir) = mpart/qpart* ( upar*vedotgradb(1:ndir))
791 call cross(drift1,drift5,vedotgradbdrift)
792 vedotgradbdrift_abs = sqrt(sum(vedotgradbdrift(:)**2))
793
794 drift6(1:ndir) = mpart/qpart* ( upar*bdotgradve(1:ndir))
795 call cross(drift1,drift6,bdotgradvedrift)
796 bdotgradvedrift_abs = sqrt(sum(bdotgradvedrift(:)**2))
797
798 drift7(1:ndir) = mpart/qpart* (upart(3)*vedotgradve(1:ndir))
799 call cross(drift1,drift7,vedotgradvedrift)
800 vedotgradvedrift_abs = sqrt(sum(vedotgradvedrift(:)**2))
801
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(:))
806
807 ! Payload update
808 if (mynpayload > 0) then
809 ! current gyroradius
810 mypayload(1) = sqrt(2.0d0*mpart*upart(2)*absb)/abs(qpart*absb)
811 end if
812 if (mynpayload > 1) then
813 ! pitch angle
814 mypayload(2) = atan2(sqrt((2.0d0*upart(2)*absb)/(mpart*upart(3)**2)),vpar)
815 end if
816 if (mynpayload > 2) then
817 ! particle v_perp
818 mypayload(3) = sqrt((2.0d0*upart(2)*absb)/(mpart*upart(3)**2))
819 end if
820 if (mynpayload > 3) then
821 ! particle parallel momentum term 1
822 mypayload(4) = momentumpar1
823 end if
824 if (mynpayload > 4) then
825 ! particle parallel momentum term 2
826 mypayload(5) = momentumpar2
827 end if
828 if (mynpayload > 5) then
829 ! particle parallel momentum term 3
830 mypayload(6) = momentumpar3
831 end if
832 if (mynpayload > 6) then
833 ! particle parallel momentum term 4
834 mypayload(7) = momentumpar4
835 end if
836 if (mynpayload > 7) then
837 ! particle ExB drift
838 mypayload(8) = veabs
839 end if
840 if (mynpayload > 8) then
841 ! relativistic drift
842 mypayload(9) = gradbdrift_abs
843 end if
844 if (mynpayload > 9) then
845 ! gradB drift
846 mypayload(10) = reldrift_abs
847 end if
848 if (mynpayload > 10) then
849 ! bdotgradb drift
850 mypayload(11) = bdotgradbdrift_abs
851 end if
852 if (mynpayload > 11) then
853 ! vEdotgradb drift
854 mypayload(12) = vedotgradbdrift_abs
855 end if
856 if (mynpayload > 12) then
857 ! bdotgradvE drift
858 mypayload(13) = bdotgradvedrift_abs
859 end if
860 if (mynpayload > 13) then
861 ! vEdotgradvE drift
862 mypayload(14) = vedotgradvedrift_abs
863 end if
864
865 end subroutine gca_update_payload
866
867 function gca_get_particle_dt(partp, end_time) result(dt_p)
868 use mod_odeint
870 type(particle_ptr), intent(in) :: partp
871 double precision, intent(in) :: end_time
872 double precision :: dt_p
873
874 double precision :: tout, dt_particles_mype, dt_cfl0, dt_cfl1, dt_a
875 double precision :: dxmin, vp, a, gammap
876 double precision :: v(ndir), y(ndir+2),ytmp(ndir+2), dydt(ndir+2), v0(ndir), v1(ndir), dydt1(ndir+2)
877 double precision :: ap0, ap1, dt_cfl_ap0, dt_cfl_ap1, dt_cfl_ap
878 double precision :: dt_euler, dt_tmp
879 ! make these particle cfl conditions more restrictive if you are interpolating out of the grid
880 double precision :: cfl, uparcfl
881 double precision, parameter :: uparmin=1.0d-6*const_c
882 integer :: ipart, iipart, nout, ic^d, igrid_particle, ipe_particle, ipe
883 logical :: bc_applied
884
885 if (const_dt_particles > 0) then
886 dt_p = const_dt_particles
887 return
888 end if
889
890 cfl = particles_cfl
891 uparcfl = particles_cfl
892
893 igrid_working = partp%igrid
894 ipart_working = partp%self%index
895 dt_tmp = (end_time - partp%self%time)
896 if(dt_tmp .le. 0.0d0) dt_tmp = smalldouble
897 ! make sure we step only one cell at a time, first check CFL at current location
898 ! then we make an Euler step to the new location and check the new CFL
899 ! we simply take the minimum of the two timesteps.
900 ! added safety factor cfl:
901 dxmin = min({rnode(rpdx^d_,igrid_working)},bigdouble)*cfl
902 ! initial solution vector:
903 y(1:ndir) = partp%self%x(1:ndir) ! position of guiding center
904 y(ndir+1) = partp%self%u(1) ! parallel momentum component (gamma v||)
905 y(ndir+2) = partp%self%u(2) ! conserved magnetic moment
906 ytmp=y
907 !y(ndir+3) = partp%self%u(3) ! Lorentz factor of guiding centre
908
909 call derivs_gca(partp%self%time,y,dydt)
910 v0(1:ndir) = dydt(1:ndir)
911 ap0 = dydt(ndir+1)
912
913 ! guiding center velocity:
914 v(1:ndir) = abs(dydt(1:ndir))
915 vp = sqrt(sum(v(:)**2))
916
917 dt_cfl0 = dxmin / max(vp, smalldouble)
918 dt_cfl_ap0 = uparcfl * abs(max(abs(y(ndir+1)),uparmin) / max(abs(ap0), smalldouble))
919 !dt_cfl_ap0 = min(dt_cfl_ap0, uparcfl * sqrt(abs(unit_length*dxmin/(ap0+smalldouble))) )
920
921 ! make an Euler step with the proposed timestep:
922 ! new solution vector:
923 dt_euler = min(dt_tmp,dt_cfl0,dt_cfl_ap0)
924 y(1:ndir+2) = y(1:ndir+2) + dt_euler * dydt(1:ndir+2)
925
926 partp%self%x(1:ndir) = y(1:ndir) ! position of guiding center
927 partp%self%u(1) = y(ndir+1) ! parallel momentum component (gamma v||)
928 partp%self%u(2) = y(ndir+2) ! conserved magnetic moment
929
930 ! first check if the particle is outside the physical domain or in the ghost cells
932 y(1:ndir+2) = ytmp(1:ndir+2)
933 end if
934
935 call derivs_gca_rk(partp%self%time+dt_euler,y,dydt)
936 !call derivs_gca(partp%self%time+dt_euler,y,dydt)
937
938 v1(1:ndir) = dydt(1:ndir)
939 ap1 = dydt(ndir+1)
940
941 ! guiding center velocity:
942 v(1:ndir) = abs(dydt(1:ndir))
943 vp = sqrt(sum(v(:)**2))
944
945 dt_cfl1 = dxmin / max(vp, smalldouble)
946 dt_cfl_ap1 = uparcfl * abs(max(abs(y(ndir+1)),uparmin) / max(abs(ap1), smalldouble))
947 !dt_cfl_ap1 = min(dt_cfl_ap1, uparcfl * sqrt(abs(unit_length*dxmin/(ap1+smalldouble))) )
948
949 dt_tmp = min(dt_euler, dt_cfl1, dt_cfl_ap1)
950
951 partp%self%x(1:ndir) = ytmp(1:ndir) ! position of guiding center
952 partp%self%u(1) = ytmp(ndir+1) ! parallel momentum component (gamma v||)
953 partp%self%u(2) = ytmp(ndir+2) ! conserved magnetic moment
954 !dt_tmp = min(dt_cfl1, dt_cfl_ap1)
955
956 ! time step due to parallel acceleration:
957 ! The standard thing, dt=sqrt(dx/a) where we compute a from d(gamma v||)/dt and d(gamma)/dt
958 ! dt_ap = sqrt(abs(dxmin*unit_length*y(ndir+3)/( dydt(ndir+1) - y(ndir+1)/y(ndir+3)*dydt(ndir+3) ) ) )
959 ! vp = sqrt(sum(v(1:ndir)**))
960 ! gammap = sqrt(1.0d0/(1.0d0-(vp/const_c)**2))
961 ! ap = const_c**2/vp*gammap**(-3)*dydt(ndir+3)
962 ! dt_ap = sqrt(dxmin*unit_length/ap)
963
964 !dt_a = bigdouble
965 !if (dt_euler .gt. smalldouble) then
966 ! a = sqrt(sum((v1(1:ndir)-v0(1:ndir))**2))/dt_euler
967 ! dt_a = min(sqrt(dxmin/a),bigdouble)
968 !end if
969
970 !dt_p = min(dt_tmp , dt_a)
971 dt_p = dt_tmp
972
973 ! Make sure we do not advance beyond end_time
974 call limit_dt_endtime(end_time - partp%self%time, dt_p)
975
976 end function gca_get_particle_dt
977
978end 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, 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 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.
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.
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