MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_particle_lorentz.t
Go to the documentation of this file.
1 !> Particle mover with Newtonian/relativistic Boris scheme for Lorentz dynamics
2 !> By Jannis Teunissen, Bart Ripperda, Oliver Porth, and Fabio Bacchini (2016-2020)
5 
6  private
7 
8  public :: lorentz_init
10  integer, parameter :: Boris=1, vay=2, hc=3, lm=4
11 
12  ! Variables
13  public :: bp, ep, vp, jp, rhop
14 
15 contains
16 
17  subroutine lorentz_init()
19  integer :: idir, nwx
20 
21  if(physics_type/='mhd') call mpistop("Lorentz particles need magnetic field!")
22  if(ndir/=3) call mpistop("Lorentz particles need ndir=3!")
23 
24  ! The first 6 gridvars are always B and E
25  ngridvars = ndir*2
26  nwx=0
27  ! density
28  if(particles_etah>0) nwx = 1
29 
30  allocate(bp(ndir))
31  do idir = 1, ndir
32  nwx = nwx + 1
33  bp(idir) = nwx
34  end do
35  allocate(ep(ndir)) ! electric field
36  do idir = 1, ndir
37  nwx = nwx + 1
38  ep(idir) = nwx
39  end do
40  allocate(vp(ndir)) ! fluid velocity
41  do idir = 1, ndir
42  nwx = nwx + 1
43  vp(idir) = nwx
44  end do
45  allocate(jp(ndir)) ! current
46  do idir = 1, ndir
47  nwx = nwx + 1
48  jp(idir) = nwx
49  end do
50  nwx = nwx + 1 ! density
51  rhop = nwx
52 
53  ngridvars=nwx
54 
56 
57  if (associated(particles_define_additional_gridvars)) then
59  end if
60 
61  select case(integrator_type_particles)
62  case('Boris','boris')
63  integrator = boris
64  case('Vay','vay')
65  integrator = vay
66  case('HC','hc','higueracary')
67  integrator = hc
68  case('LM','lm','lapentamarkidis')
69  integrator = lm
70  case default
71  integrator = boris
72  end select
73 
75 
76  end subroutine lorentz_init
77 
79 
82 
83  integer :: n, idir, igrid, ipe_particle, nparticles_local
84  double precision :: lfac
85  double precision :: x(3, num_particles)
86  double precision :: v(3, num_particles)
87  double precision :: q(num_particles)
88  double precision :: m(num_particles)
89  double precision :: rrd(num_particles,ndir)
90  double precision :: defpayload(ndefpayload)
91  double precision :: usrpayload(nusrpayload)
92  logical :: follow(num_particles), check
93 
94  follow = .false.
95  x = 0.0d0
96 
97  if (mype==0) then
98  if (.not. associated(usr_create_particles)) then
99  ! Randomly distributed
100  do idir=1,ndir
101  do n = 1, num_particles
102  rrd(n,idir) = rng%unif_01()
103  end do
104  end do
105  do n=1, num_particles
106  {^d&x(^d,n) = xprobmin^d + rrd(n+1,^d) * (xprobmax^d - xprobmin^d)\}
107  end do
108  else
109  call usr_create_particles(num_particles, x, v, q, m, follow)
110  end if
111  end if
112 
113  call mpi_bcast(x,3*num_particles,mpi_double_precision,0,icomm,ierrmpi)
114  call mpi_bcast(v,3*num_particles,mpi_double_precision,0,icomm,ierrmpi)
115  call mpi_bcast(q,num_particles,mpi_double_precision,0,icomm,ierrmpi)
116  call mpi_bcast(m,num_particles,mpi_double_precision,0,icomm,ierrmpi)
117  call mpi_bcast(follow,num_particles,mpi_logical,0,icomm,ierrmpi)
118 
119  nparticles_local = 0
120 
121  ! Find ipe and igrid responsible for particle
122  do n = 1, num_particles
123  call find_particle_ipe(x(:, n),igrid,ipe_particle)
124 
125  particle(n)%igrid = igrid
126  particle(n)%ipe = ipe_particle
127 
128  if (ipe_particle == mype) then
129  check = .true.
130 
131  ! Check for user-defined modifications or rejection conditions
132  if (associated(usr_check_particle)) call usr_check_particle(igrid, x(:,n), v(:,n), q(n), m(n), follow(n), check)
133  if (check) then
135  else
136  cycle
137  end if
138 
139  nparticles_local = nparticles_local + 1
140 
141  call get_lfac_from_velocity(v(:,n), lfac)
142 
143  allocate(particle(n)%self)
144  particle(n)%self%x = x(:,n)
145  particle(n)%self%u = v(:,n) * lfac
146  particle(n)%self%q = q(n)
147  particle(n)%self%m = m(n)
148  particle(n)%self%follow = follow(n)
149  particle(n)%self%index = n
150  particle(n)%self%time = global_time
151  particle(n)%self%dt = 0.0d0
152 
153  ! initialise payloads for Lorentz module
154  allocate(particle(n)%payload(npayload))
155  call lorentz_update_payload(igrid,x(:,n),v(:,n)*lfac,q(n),m(n),defpayload,ndefpayload,0.d0)
156  particle(n)%payload(1:ndefpayload) = defpayload
157  if (associated(usr_update_payload)) then
158  call usr_update_payload(igrid,x(:,n),v(:,n)*lfac,q(n),m(n),usrpayload,nusrpayload,0.d0)
159  particle(n)%payload(ndefpayload+1:npayload) = usrpayload
160  end if
161  end if
162  end do
163 
164  call mpi_allreduce(nparticles_local,nparticles,1,mpi_integer,mpi_sum,icomm,ierrmpi)
165 
166  end subroutine lorentz_create_particles
167 
171  use mod_geometry
172 
173  integer :: igrid, iigrid, idir, idim
174  double precision, dimension(ixG^T,1:ndir) :: beta
175  double precision, dimension(ixG^T,1:nw) :: w,wold
176  double precision :: current(ixG^T,7-2*ndir:3)
177  integer :: idirmin
178  double precision, dimension(ixG^T,1:ndir) :: vE, bhat
179  double precision, dimension(ixG^T) :: kappa, kappa_B, absB, tmp
180 
181  ! Fill electromagnetic quantities
182  call fill_gridvars_default()
183 
184  end subroutine lorentz_fill_gridvars
185 
186  !> Relativistic particle integrator
187  subroutine lorentz_integrate_particles(end_time)
189  use mod_geometry
191  double precision, intent(in) :: end_time
192  double precision :: defpayload(ndefpayload)
193  double precision :: usrpayload(nusrpayload)
194  integer :: ipart, iipart, igrid
195  double precision :: lfac, q, m, dt_p
196  double precision :: xp(ndir), xpm(ndir), xpc(ndir), xpcm(ndir)
197  double precision :: up(ndir), upc(ndir), tp
198  double precision, dimension(ndir) :: b, e, bc, ec, vfluid, current
199  double precision :: rho, rhoold, td
200 
201  do iipart=1,nparticles_active_on_mype
202  ipart = particles_active_on_mype(iipart);
203  q = particle(ipart)%self%q
204  m = particle(ipart)%self%m
205  igrid = particle(ipart)%igrid
206  xp = particle(ipart)%self%x
207  up = particle(ipart)%self%u
208  tp = particle(ipart)%self%time
209 
210  dt_p = lorentz_get_particle_dt(particle(ipart), end_time)
211  particle(ipart)%self%dt = dt_p
212 
213  ! Push particle over first half time step
214  ! all pushes done in Cartesian coordinates
215  call partcoord_to_cartesian(xp,xpc)
216  call partvec_to_cartesian(xp,up,upc)
217  call get_lfac(upc,lfac)
218  xpcm = xpc + 0.5d0 * dt_p * upc/lfac
219  call partcoord_from_cartesian(xpm,xpcm)
220  ! Fix xp if the 0,2*pi boundary was crossed in cylindrical/spherical coords
221  call fix_phi_crossing(xpm,igrid)
222 
223  ! Get E, B at n+1/2 position
224  call get_bfield(igrid, xpm, tp+dt_p/2.d0, b)
225  call get_efield(igrid, xpm, tp+dt_p/2.d0, e)
226 
227 ! call get_vec(vp, igrid, xpm, tp+dt_p/2.d0, vfluid)
228 ! call get_vec(jp, igrid, xpm, tp+dt_p/2.d0, current)
229 ! select case (coordinate)
230 ! case (Cartesian,Cartesian_stretched,spherical)
231 ! e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
232 ! e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
233 ! e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)
234 ! case (cylindrical)
235 ! e(r_) = -vfluid(phi_)*b(z_)+vfluid(z_)*b(phi_) + particles_eta*current(r_)
236 ! e(phi_) = vfluid(r_)*b(z_)-vfluid(z_)*b(r_) + particles_eta*current(phi_)
237 ! e(z_) = -vfluid(r_)*b(phi_)+vfluid(phi_)*b(r_) + particles_eta*current(z_)
238 ! end select
239 ! if (particles_etah > zero) then
240 ! call interpolate_var(igrid,ixG^LL,ixM^LL,ps(igrid)%w(ixG^T,1),ps(igrid)%x,xpm,rho)
241 ! if (time_advance) then
242 ! td = (tp+dt_p/2.d0 - global_time) / dt
243 ! call interpolate_var(igrid,ixG^LL,ixM^LL,pso(igrid)%w(ixG^T,1),ps(igrid)%x,xpm,rhoold)
244 ! rho = rhoold * (1.d0-td) + rho * td
245 ! end if
246 ! select case (coordinate)
247 ! case (Cartesian,Cartesian_stretched,spherical)
248 ! e(1) = e(1) + particles_etah/rho * (current(2)*b(3) - current(3)*b(2))
249 ! e(2) = e(2) + particles_etah/rho * (-current(1)*b(3) + current(3)*b(1))
250 ! e(3) = e(3) + particles_etah/rho * (current(1)*b(2) - current(2)*b(1))
251 ! case (cylindrical)
252 ! e(r_) = e(r_) + particles_etah/rho * (current(phi_)*b(z_) - current(z_)*b(phi_))
253 ! e(phi_) = e(phi_) + particles_etah/rho * (-current(r_)*b(z_) + current(z_)*b(r_))
254 ! e(z_) = e(z_) + particles_etah/rho * (current(r_)*b(phi_) - current(phi_)*b(r_))
255 ! end select
256 ! end if
257 
258  ! Convert fields to Cartesian frame
259  call partvec_to_cartesian(xpm,b,bc)
260  call partvec_to_cartesian(xpm,e,ec)
261 
262  ! 'Kick' particle (update velocity) based on the chosen integrator
263  call lorentz_kick(upc,ec,bc,q,m,dt_p)
264 
265  ! Push particle over second half time step
266  ! all pushes done in Cartesian coordinates
267  call get_lfac(upc,lfac)
268  xpc = xpcm + 0.5d0 * dt_p * upc/lfac
269  call partcoord_from_cartesian(xp,xpc)
270  ! Fix xp if the 0,2*pi boundary was crossed in cylindrical/spherical coords
271  call fix_phi_crossing(xp,igrid)
272  call partvec_from_cartesian(xp,up,upc)
273 
274  ! Store updated x,u
275  particle(ipart)%self%x = xp
276  particle(ipart)%self%u = up
277 
278  ! Time update
279  tp = tp + dt_p
280  particle(ipart)%self%time = tp
281 
282  ! Update payload
283  call lorentz_update_payload(igrid,xp,up,q,m,defpayload,ndefpayload,tp)
284  particle(ipart)%payload(1:ndefpayload) = defpayload
285  if (associated(usr_update_payload)) then
286  call usr_update_payload(igrid,xp,up,q,m,usrpayload,nusrpayload,tp)
287  particle(ipart)%payload(ndefpayload+1:npayload) = usrpayload
288  end if
289 
290  end do ! ipart loop
291 
292  end subroutine lorentz_integrate_particles
293 
294  !> Momentum update subroutine for full Lorentz dynamics
295  subroutine lorentz_kick(upart,e,b,q,m,dtp)
297  use mod_geometry
298  double precision, intent(in) :: e(ndir), b(ndir), q, m, dtp
299  double precision, intent(inout) :: upart(ndir)
300  double precision :: lfac, cosphi, sinphi, phi1, phi2, r, re, sigma, td
301  double precision, dimension(ndir) :: emom, uprime, tau, s, tmp, uplus, xcart1, xcart2, ucart2, radmom
302  double precision, dimension(ndir) :: upartk, vbar, Fk, C1, C2, dupartk
303  double precision :: abserr, tol, lfack, J11, J12, J13, J21, J22, J23, J31, J32, J33
304  double precision :: iJ11, iJ12, iJ13, iJ21, iJ22, iJ23, iJ31, iJ32, iJ33, Det
305  integer :: nk, nkmax
306 
307  ! Perform momentum update based on the chosen integrator
308  select case(integrator)
309 
310  ! Boris integrator (works in Cartesian and cylindrical)
311  case(boris)
312  ! Momentum update
313  emom = q * e * dtp /(2.0d0 * m)
314  uprime = upart + emom
315  !!!!!! TODO: Adjust and document losses
316 ! if (losses) then
317 ! call get_lfac(particle(ipart)%self%u,lfac)
318 ! re = abs(q)**2 / (m * const_c**2)
319 ! call cross(upart,b,tmp)
320 ! radmom = - third * re**2 * lfac &
321 ! * ( sum((e(:)+tmp(:)/lfac)**2) &
322 ! - (sum(e(:)*upart(:))/lfac)**2 ) &
323 ! * particle(ipart)%self%u / m / const_c * dt_p
324 ! uprime = uprime + radmom
325 ! end if
326 
327  call get_lfac(uprime,lfac)
328  tau = q * b * dtp / (2.0d0 * lfac * m)
329  s = 2.0d0 * tau / (1.0d0+sum(tau(:)**2))
330 
331  call cross(uprime,tau,tmp)
332  call cross(uprime+tmp,s,tmp)
333  uplus = uprime + tmp
334 
335  upart = uplus + emom
336  !!!!!! TODO: Adjust and document losses
337 ! if(losses) then
338 ! call cross(uplus,b,tmp)
339 ! radmom = - third * re**2 * lfac &
340 ! * ( sum((e(:)+tmp(:)/lfac)**2) &
341 ! - (sum(e(:)*uplus(:))/lfac)**2 ) &
342 ! * uplus / m / const_c * dt_p
343 ! upart = upart + radmom
344 ! end if
345 
346  ! Vay integrator
347  case(vay)
348  call get_lfac(upart,lfac)
349  emom = q * e * dtp /(2.0d0 * m)
350  tau = q * b * dtp / (2.0d0 * m)
351 
352  call cross(upart,tau,tmp)
353  uprime = upart + 2.d0*emom + tmp/lfac
354 
355  call get_lfac(uprime,lfac)
356  sigma = lfac**2 - sum(tau(:)*tau(:))
357  lfac = sqrt((sigma + sqrt(sigma**2 &
358  + 4.d0 * (sum(tau(:)*tau(:)) + sum(uprime(:)*tau(:)/c_norm)**2))) / 2.d0)
359 
360  call cross(uprime,tau,tmp)
361  upart = (uprime + sum(uprime(:)*tau(:))*tau/lfac**2 + tmp/lfac) / (1.d0+sum(tau(:)*tau(:))/lfac**2)
362 
363  ! Higuera-Cary integrator
364  case(hc)
365  call get_lfac(upart,lfac)
366  emom = q * e * dtp /(2.0d0 * m)
367  tau = q * b * dtp / (2.0d0 * m)
368  uprime = upart + emom
369 
370  call get_lfac(uprime,lfac)
371  sigma = lfac**2 - sum(tau(:)*tau(:))
372  lfac = sqrt((sigma + sqrt(sigma**2 &
373  + 4.d0 * (sum(tau(:)*tau(:)) + sum(uprime(:)*tau(:)/c_norm)**2))) / 2.d0)
374 
375  call cross(uprime,tau,tmp)
376  upart = (uprime + sum(uprime(:)*tau(:))*tau/lfac**2 + tmp/lfac) / (1.d0+sum(tau(:)*tau(:))/lfac**2) &
377  + emom + tmp/lfac
378 
379  ! Lapenta-Markidis integrator
380  case(lm)
381  ! Initialise iteration quantities
382  call get_lfac(upart,lfac)
383  upartk = upart
384  tau(:) = b(:)
385 
386  ! START OF THE NONLINEAR CYCLE
387  abserr = 1.d0
388  tol=1.d-14
389  nkmax=10
390  nk=0
391  do while(abserr > tol .and. nk < nkmax)
392 
393  nk=nk+1
394 
395  call get_lfac(upartk,lfack)
396  vbar = (upart + upartk) / (lfac + lfack)
397  call cross(vbar,tau,tmp)
398 
399  ! Compute residual vector
400  fk = upartk - upart - q*dtp/m * (e + tmp)
401 
402  ! Compute auxiliary coefficients
403  c1 = (lfack + lfac - upartk(1:ndim) / lfack / c_norm**2 * (upartk + upart)) / (lfack + lfac)**2
404  c2 = - upartk / lfack / c_norm**2 / (lfack + lfac)**2
405 
406  ! Compute Jacobian
407  j11 = 1. - q*dtp/m * (c2(1) * (upartk(2) + upart(2)) * tau(3) - c2(1) * (upartk(3) + upart(3)) * tau(2))
408  j12 = - q*dtp/m * (c1(2) * tau(3) - c2(2) * (upartk(3) + upart(3)) * tau(2))
409  j13 = - q*dtp/m * (c2(3) * (upartk(2) + upart(2)) * tau(3) - c1(3) * tau(2))
410  j21 = - q*dtp/m * (- c1(1) * tau(3) + c2(1) * (upartk(3) + upart(3)) * tau(1))
411  j22 = 1. - q*dtp/m * (- c2(2) * (upartk(1) + upart(1)) * tau(3) + c2(2) * (upartk(3) + upart(3)) * tau(1))
412  j23 = - q*dtp/m * (- c2(3) * (upartk(1) + upart(1)) * tau(3) + c1(3) * tau(1))
413  j31 = - q*dtp/m * (c1(1) * tau(2) - c2(1) * (upartk(2) + upart(2)) * tau(1))
414  j32 = - q*dtp/m * (c2(2) * (upartk(1) + upart(1)) * tau(2) - c1(2) * tau(1))
415  j33 = 1. - q*dtp/m * (c2(3) * (upartk(1) + upart(1)) * tau(2) - c2(3) * (upartk(2) + upart(2)) * tau(1))
416 
417  ! Compute inverse Jacobian
418  det = j11*j22*j33 + j21*j32*j13 + j31*j12*j23 - j11*j32*j23 - j31*j22*j13 - j21*j12*j33
419  ij11 = (j22*j33 - j23*j32) / det
420  ij12 = (j13*j32 - j12*j33) / det
421  ij13 = (j12*j23 - j13*j22) / det
422  ij21 = (j23*j31 - j21*j33) / det
423  ij22 = (j11*j33 - j13*j31) / det
424  ij23 = (j13*j21 - j11*j23) / det
425  ij31 = (j21*j32 - j22*j31) / det
426  ij32 = (j12*j31 - j11*j32) / det
427  ij33 = (j11*j22 - j12*j21) / det
428 
429  ! Compute new upartk = upartk - J^(-1) * F(upartk)
430  dupartk(1) = - (ij11 * fk(1) + ij12 * fk(2) + ij13 * fk(3))
431  dupartk(2) = - (ij21 * fk(1) + ij22 * fk(2) + ij23 * fk(3))
432  dupartk(3) = - (ij31 * fk(1) + ij32 * fk(2) + ij33 * fk(3))
433 
434  ! Check convergence
435  upartk=upartk+dupartk
436  abserr=sqrt(sum(dupartk(:)*dupartk(:)))
437 
438  end do
439  ! END OF THE NONLINEAR CYCLE
440 
441  ! Update velocity
442  upart = upartk
443 
444  end select
445 
446  end subroutine lorentz_kick
447 
448  !> Update payload subroutine
449  subroutine lorentz_update_payload(igrid,xpart,upart,qpart,mpart,mypayload,mynpayload,particle_time)
451  use mod_geometry
452  integer, intent(in) :: igrid,mynpayload
453  double precision, intent(in) :: xpart(1:ndir),upart(1:ndir),qpart,mpart,particle_time
454  double precision, intent(out) :: mypayload(mynpayload)
455  double precision :: b(3), e(3), tmp(3), lfac, vfluid(3), current(3), rho, rhoold, td
456 
457  call get_bfield(igrid, xpart, particle_time, b)
458  call get_efield(igrid, xpart, particle_time, e)
459 
460 ! call get_vec(bp, igrid, xpart,particle_time,b)
461 ! call get_vec(vp, igrid, xpart,particle_time,vfluid)
462 ! call get_vec(jp, igrid, xpart,particle_time,current)
463 ! select case (coordinate)
464 ! case (Cartesian,Cartesian_stretched,spherical)
465 ! e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
466 ! e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
467 ! e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)
468 ! case (cylindrical)
469 ! e(r_) = -vfluid(phi_)*b(z_)+vfluid(z_)*b(phi_) + particles_eta*current(r_)
470 ! e(phi_) = vfluid(r_)*b(z_)-vfluid(z_)*b(r_) + particles_eta*current(phi_)
471 ! e(z_) = -vfluid(r_)*b(phi_)+vfluid(phi_)*b(r_) + particles_eta*current(z_)
472 ! end select
473 ! if (particles_etah > zero) then
474 ! call interpolate_var(igrid,ixG^LL,ixM^LL,ps(igrid)%w(ixG^T,1),ps(igrid)%x,xpart,rho)
475 ! if (time_advance) then
476 ! td = (particle_time - global_time) / dt
477 ! call interpolate_var(igrid,ixG^LL,ixM^LL,pso(igrid)%w(ixG^T,1),ps(igrid)%x,xpart,rhoold)
478 ! rho = rhoold * (1.d0-td) + rho * td
479 ! end if
480 ! select case (coordinate)
481 ! case (Cartesian,Cartesian_stretched,spherical)
482 ! e(1) = e(1) + particles_etah/rho * (current(2)*b(3) - current(3)*b(2))
483 ! e(2) = e(2) + particles_etah/rho * (-current(1)*b(3) + current(3)*b(1))
484 ! e(3) = e(3) + particles_etah/rho * (current(1)*b(2) - current(2)*b(1))
485 ! case (cylindrical)
486 ! e(r_) = e(r_) + particles_etah/rho * (current(phi_)*b(z_) - current(z_)*b(phi_))
487 ! e(phi_) = e(phi_) + particles_etah/rho * (-current(r_)*b(z_) + current(z_)*b(r_))
488 ! e(z_) = e(z_) + particles_etah/rho * (current(r_)*b(phi_) - current(phi_)*b(r_))
489 ! end select
490 ! end if
491 
492  ! Payload update
493  ! Lorentz factor
494  if (mynpayload > 0) then
495  call get_lfac(upart,lfac)
496  mypayload(1) = lfac
497  end if
498 
499  ! current gyroradius
500  if (mynpayload > 1) then
501  call cross(upart,b,tmp)
502  tmp = tmp / sqrt(sum(b(:)**2))
503  mypayload(2) = mpart/abs(qpart)*sqrt(sum(tmp(:)**2)) / sqrt(sum(b(:)**2))
504  end if
505 
506  ! magnetic moment
507  if (mynpayload > 2) then
508  mypayload(3) = mpart*sum(tmp(:)**2)/(2.0d0*sqrt(sum(b(:)**2)))
509  end if
510 
511  ! e.b
512  if (mynpayload > 3) then
513  mypayload(4) = sum(e(:)*b(:))
514  end if
515 
516  end subroutine lorentz_update_payload
517 
518  function lorentz_get_particle_dt(partp, end_time) result(dt_p)
520  use mod_geometry
521  type(particle_ptr), intent(in) :: partp
522  double precision, intent(in) :: end_time
523  double precision :: dt_p
524  integer :: ipart, iipart, nout
525  double precision,dimension(ndir) :: b,v
526  double precision :: lfac,absb,dt_cfl
527  double precision :: tout
528  double precision, parameter :: cfl = 0.5d0
529 
530  if (const_dt_particles > 0) then
531  dt_p = const_dt_particles
532  return
533  end if
534 
535  call get_vec(bp, partp%igrid,partp%self%x,partp%self%time,b)
536  absb = sqrt(sum(b(:)**2))
537  call get_lfac(partp%self%u,lfac)
538 
539  ! CFL timestep
540  ! make sure we step only one cell at a time:
541  v(:) = abs(partp%self%u(:) / lfac)
542 
543 ! ! convert to angular velocity:
544 ! if(coordinate ==cylindrical.and.phi_>0) then
545 ! v(phi_) = abs(v(phi_)/partp%self%x(r_))
546 ! end if
547 
548  dt_cfl = min(bigdouble, &
549  {rnode(rpdx^d_,partp%igrid)/max(v(^d), smalldouble)})
550 
551  if(coordinate ==cylindrical.and.phi_>0) then
552  ! phi-momentum leads to radial velocity:
553  if(phi_ .gt. ndim) dt_cfl = min(dt_cfl, &
554  sqrt(rnode(rpdx1_,partp%igrid)/partp%self%x(r_)) &
555  / v(phi_))
556  ! limit the delta phi of the orbit (just for aesthetic reasons):
557  dt_cfl = min(dt_cfl,0.1d0/max(v(phi_), smalldouble))
558  ! take some care at the axis:
559  dt_cfl = min(dt_cfl,(partp%self%x(r_)+smalldouble)/max(v(r_), smalldouble))
560  end if
561 
562  dt_cfl = dt_cfl * cfl
563 
564  ! bound by gyro-rotation:
565  dt_p = abs( dtheta * partp%self%m * lfac &
566  / (partp%self%q * absb) )
567 
568  dt_p = min(dt_p, dt_cfl)
569 
570  ! Make sure we don't advance beyond end_time
571  call limit_dt_endtime(end_time - partp%self%time, dt_p)
572 
573  end function lorentz_get_particle_dt
574 
575 end module mod_particle_lorentz
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
integer coordinate
Definition: mod_geometry.t:7
integer, parameter cylindrical
Definition: mod_geometry.t:10
This module contains definitions of global parameters and variables and some generic functions/subrou...
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 r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
Module with shared functionality for all the particle movers.
subroutine partcoord_from_cartesian(xp, xpcart)
Convert to noncartesian coordinates.
pure subroutine limit_dt_endtime(t_left, dt_p)
integer num_particles
Number of particles.
subroutine fix_phi_crossing(xp, igrid)
Fix particle position when crossing the 0,2pi boundary in noncartesian coordinates.
double precision const_dt_particles
If positive, a constant time step for the particles.
integer ngridvars
Number of variables for grid field.
subroutine get_efield(igrid, x, tloc, e)
double precision particles_etah
integer nusrpayload
Number of user-defined payload variables for a particle.
pure subroutine get_lfac(u, lfac)
Get Lorentz factor from relativistic momentum.
subroutine find_particle_ipe(x, igrid_particle, ipe_particle)
double precision dtheta
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 partvec_from_cartesian(xp, up, upcart)
Convert vector from Cartesian coordinates.
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.
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
subroutine partcoord_to_cartesian(xp, xpcart)
Convert position to Cartesian coordinates.
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
subroutine partvec_to_cartesian(xp, up, upcart)
Convert vector to Cartesian coordinates.
integer ndefpayload
Number of default payload variables for a particle.
integer, dimension(:), allocatable jp
Variable index for current.
integer rhop
Variable index for density.
Particle mover with Newtonian/relativistic Boris scheme for Lorentz dynamics By Jannis Teunissen,...
double precision function lorentz_get_particle_dt(partp, end_time)
subroutine lorentz_update_payload(igrid, xpart, upart, qpart, mpart, mypayload, mynpayload, particle_time)
Update payload subroutine.
subroutine, public lorentz_init()
subroutine lorentz_integrate_particles(end_time)
Relativistic particle integrator.
subroutine lorentz_kick(upart, e, b, q, m, dtp)
Momentum update subroutine for full Lorentz dynamics.
subroutine, public lorentz_create_particles()
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