36 if (physics_type/=
'mhd')
call mpistop(
"GCA particles need magnetic field!")
37 if (
ndir/=3)
call mpistop(
"GCA particles need ndir=3!")
102 double precision :: b(
ndir), u(
ndir), magmom
103 double precision :: bnorm, lfac, vnorm, vperp, vpar
104 integer :: igrid, ipe_particle
105 integer :: n, idir, nparticles_local
120 rrd(n,idir) =
rng%unif_01()
124 {^
d&x(^
d,n) = xprobmin^
d + rrd(n+1,^
d) * (xprobmax^
d - xprobmin^
d)\}
146 if(ipe_particle ==
mype)
then
157 nparticles_local = nparticles_local + 1
173 vnorm = norm2(v(:, n))
174 vpar = sum(v(:, n) * b/bnorm)
175 vperp = sqrt(vnorm**2 - vpar**2)
183 magmom = m(n) * (vperp * lfac)**2 / (2.0d0 * bnorm)
191 call gca_update_payload(igrid,ps(igrid)%w,pso(igrid)%w,ps(igrid)%x,x(:,n),
particle(n)%self%u(:),q(n),m(n),defpayload,
ndefpayload,0.d0)
194 call usr_update_payload(igrid,ps(igrid)%w,pso(igrid)%w,ps(igrid)%x,x(:,n),
particle(n)%self%u(:),q(n),m(n),usrpayload,
nusrpayload,0.d0)
209 integer :: igrid, iigrid, idir, idim
210 double precision,
dimension(ixG^T,1:ndir) :: beta
211 double precision,
dimension(ixG^T,1:nw) :: w,wold
212 double precision :: current(ixG^T,7-2*ndir:3)
214 double precision,
dimension(ixG^T,1:ndir) :: vE, bhat
215 double precision,
dimension(ixG^T) :: kappa, kappa_B, absB, tmp
219 do iigrid=1,igridstail; igrid=igrids(iigrid);
220 w(ixg^t,1:nw) = ps(igrid)%w(ixg^t,1:nw)
221 call phys_to_primitive(ixg^
ll,ixg^
ll,w,ps(igrid)%x)
223 gridvars(igrid)%w(ixg^t,
vp(:)) = w(ixg^t,iw_mom(:))
226 absb(ixg^t) = sqrt(sum(gridvars(igrid)%w(ixg^t,
bp(:))**2,dim=
ndim+1))
227 ve(ixg^t,1) = gridvars(igrid)%w(ixg^t,
ep(2)) * gridvars(igrid)%w(ixg^t,
bp(3)) &
228 - gridvars(igrid)%w(ixg^t,
ep(3)) * gridvars(igrid)%w(ixg^t,
bp(2))
229 ve(ixg^t,2) = gridvars(igrid)%w(ixg^t,
ep(3)) * gridvars(igrid)%w(ixg^t,
bp(1)) &
230 - gridvars(igrid)%w(ixg^t,
ep(1)) * gridvars(igrid)%w(ixg^t,
bp(3))
231 ve(ixg^t,3) = gridvars(igrid)%w(ixg^t,
ep(1)) * gridvars(igrid)%w(ixg^t,
bp(2)) &
232 - gridvars(igrid)%w(ixg^t,
ep(2)) * gridvars(igrid)%w(ixg^t,
bp(1))
234 where (absb(ixg^t) .gt. 0.d0)
235 ve(ixg^t,idir) = ve(ixg^t,idir) / absb(ixg^t)**2
237 ve(ixg^t,idir) = 0.d0
241 call mpistop(
"GCA FILL GRIDVARS: vE>c! ABORTING...")
243 if (any(ve .ne. ve))
then
244 call mpistop(
"GCA FILL GRIDVARS: NaNs IN vE! ABORTING...")
248 kappa(ixg^t) = 1.d0/sqrt(1.0d0 - sum(ve(ixg^t,:)**2,dim=
ndim+1)/
c_norm**2)
252 kappa_b(ixg^t) = absb(ixg^t) / kappa(ixg^t)
254 if (any(kappa_b .ne. kappa_b))
then
255 call mpistop(
"GCA FILL GRIDVARS: NaNs IN kappa_B! ABORTING...")
260 gridvars(igrid)%w(ixg^t,
grad_kappa_b(idim)) = tmp(ixg^t)
265 where (absb(ixg^t) .gt. 0.d0)
266 bhat(ixg^t,idir) = gridvars(igrid)%w(ixg^t,
bp(idir)) / absb(ixg^t)
268 bhat(ixg^t,idir) = 0.d0
271 if (any(bhat .ne. bhat))
then
272 call mpistop(
"GCA FILL GRIDVARS: NaNs IN bhat! ABORTING...")
278 call gradient(bhat(ixg^t,idir),ixg^
ll,ixg^
ll^lsub1,idim,tmp)
280 + bhat(ixg^t,idim) * tmp(ixg^t)
282 + ve(ixg^t,idim) * tmp(ixg^t)
283 call gradient(ve(ixg^t,idir),ixg^
ll,ixg^
ll^lsub1,idim,tmp)
285 + bhat(ixg^t,idim) * tmp(ixg^t)
287 + ve(ixg^t,idim) * tmp(ixg^t)
293 w(ixg^t,1:nw) = pso(igrid)%w(ixg^t,1:nw)
294 call phys_to_primitive(ixg^
ll,ixg^
ll,w,ps(igrid)%x)
295 gridvars(igrid)%wold(ixg^t,
vp(:)) = w(ixg^t,iw_mom(:))
298 absb(ixg^t) = sqrt(sum(gridvars(igrid)%wold(ixg^t,
bp(:))**2,dim=
ndim+1))
299 ve(ixg^t,1) = gridvars(igrid)%wold(ixg^t,
ep(2)) * gridvars(igrid)%wold(ixg^t,
bp(3)) &
300 - gridvars(igrid)%wold(ixg^t,
ep(3)) * gridvars(igrid)%wold(ixg^t,
bp(2))
301 ve(ixg^t,2) = gridvars(igrid)%wold(ixg^t,
ep(3)) * gridvars(igrid)%wold(ixg^t,
bp(1)) &
302 - gridvars(igrid)%wold(ixg^t,
ep(1)) * gridvars(igrid)%wold(ixg^t,
bp(3))
303 ve(ixg^t,3) = gridvars(igrid)%wold(ixg^t,
ep(1)) * gridvars(igrid)%wold(ixg^t,
bp(2)) &
304 - gridvars(igrid)%wold(ixg^t,
ep(2)) * gridvars(igrid)%wold(ixg^t,
bp(1))
306 where (absb(ixg^t) .gt. 0.d0)
307 ve(ixg^t,idir) = ve(ixg^t,idir) / absb(ixg^t)**2
309 ve(ixg^t,idir) = 0.d0
313 call mpistop(
"GCA FILL GRIDVARS: vE>c! ABORTING...")
315 if (any(ve .ne. ve))
then
316 call mpistop(
"GCA FILL GRIDVARS: NaNs IN vE! ABORTING...")
320 kappa(ixg^t) = 1.d0/sqrt(1.0d0 - sum(ve(ixg^t,:)**2,dim=
ndim+1)/
c_norm**2)
324 kappa_b(ixg^t) = absb(ixg^t) / kappa(ixg^t)
325 if (any(kappa_b .ne. kappa_b))
then
326 call mpistop(
"GCA FILL GRIDVARS: NaNs IN kappa_B! ABORTING...")
331 gridvars(igrid)%wold(ixg^t,
grad_kappa_b(idim)) = tmp(ixg^t)
335 where (absb(ixg^t) .gt. 0.d0)
336 bhat(ixg^t,idir) = gridvars(igrid)%wold(ixg^t,
bp(idir)) / absb(ixg^t)
338 bhat(ixg^t,idir) = 0.d0
341 if (any(bhat .ne. bhat))
then
342 call mpistop(
"GCA FILL GRIDVARS: NaNs IN bhat! ABORTING...")
348 call gradient(bhat(ixg^t,idir),ixg^
ll,ixg^
ll^lsub1,idim,tmp)
350 + bhat(ixg^t,idim) * tmp(ixg^t)
352 + ve(ixg^t,idim) * tmp(ixg^t)
353 call gradient(ve(ixg^t,idir),ixg^
ll,ixg^
ll^lsub1,idim,tmp)
355 + bhat(ixg^t,idim) * tmp(ixg^t)
357 + ve(ixg^t,idim) * tmp(ixg^t)
369 double precision,
intent(in) :: end_time
371 double precision :: lfac, absS
372 double precision :: defpayload(ndefpayload)
373 double precision :: usrpayload(nusrpayload)
374 double precision :: dt_p, tloc, y(ndir+2),dydt(ndir+2),ytmp(ndir+2), euler_cfl, int_factor
375 double precision,
dimension(1:ndir) :: x, vE, e, b, bhat, x_new, vfluid, current
376 double precision,
dimension(1:ndir) :: drift1, drift2
377 double precision,
dimension(1:ndir) :: drift3, drift4, drift5, drift6, drift7
378 double precision,
dimension(1:ndir) :: bdotgradb, vEdotgradb, gradkappaB
379 double precision,
dimension(1:ndir) :: bdotgradvE, vEdotgradvE
380 double precision,
dimension(1:ndir) :: gradBdrift, reldrift, bdotgradbdrift
381 double precision,
dimension(1:ndir) :: vEdotgradbdrift, bdotgradvEdrift
382 double precision,
dimension(1:ndir) :: vEdotgradvEdrift
383 double precision :: kappa, Mr, upar, m, absb, gamma, q, mompar, vpar, vEabs
384 double precision :: gradBdrift_abs, reldrift_abs, epar
385 double precision :: bdotgradbdrift_abs, vEdotgradbdrift_abs
386 double precision :: bdotgradvEdrift_abs, vEdotgradvEdrift_abs
387 double precision :: momentumpar1, momentumpar2, momentumpar3, momentumpar4
389 double precision,
parameter :: eps=1.0
d-6
391 double precision :: h1, hmin, h_old
392 integer :: nok, nbad, ic1^D, ic2^D, ierror, nvar
393 integer :: ipart, iipart, seed, ic^D,igrid_particle, ipe_particle, ipe_working
394 logical :: int_choice
395 logical :: BC_applied
408 x(1:ndir) =
particle(ipart)%self%x(1:ndir)
412 y(1:ndir) = x(1:ndir)
413 y(ndir+1) =
particle(ipart)%self%u(1)
414 y(ndir+2) =
particle(ipart)%self%u(2)
429 y(1:ndir+2) = y(1:ndir+2) + euler_cfl * dt_p * dydt(1:ndir+2)
430 particle(ipart)%self%x(1:ndir) = y(1:ndir)
431 particle(ipart)%self%u(1) = y(ndir+1)
432 particle(ipart)%self%u(2) = y(ndir+2)
461 if (ic1^d .le.
ixglo^d-2 .or. ic2^d .ge.
ixghi^d+2)
then
467 if (ic1^d .eq.
ixglo^d-1 .or. ic2^d .eq.
ixghi^d+1)
then
476 y(1:ndir+2) = ytmp(1:ndir+2)
478 particle(ipart)%self%x(1:ndir) = ytmp(1:ndir)
479 particle(ipart)%self%u(1) = ytmp(ndir+1)
480 particle(ipart)%self%u(2) = ytmp(ndir+2)
484 h1 = dt_p/2.0d0; hmin=1.0d-9; h_old=dt_p/2.0d0
491 if (any(y .ne. y))
then
492 call mpistop(
"NaNs DETECTED IN GCA_INTEGRATE BEFORE ODEINT CALL! ABORTING...")
496 call odeint(y,nvar,tloc,tloc+dt_p,eps,h1,hmin,nok,nbad,
derivs_gca_rk,
rkqs,ierror)
498 if (ierror /= 0)
then
499 print *,
"odeint returned error code", ierror
500 print *,
"1 means hmin too small, 2 means MAXSTP exceeded"
501 print *,
"Having a problem with particle", iipart
504 if (any(y .ne. y))
then
505 call mpistop(
"NaNs DETECTED IN GCA_INTEGRATE AFTER ODEINT CALL! ABORTING...")
512 particle(ipart)%self%x(1:ndir) = y(1:ndir)
513 particle(ipart)%self%u(1) = y(ndir+1)
514 particle(ipart)%self%u(2) = y(ndir+2)
521 e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) +
particles_eta*current(1)
522 e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) +
particles_eta*current(2)
523 e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) +
particles_eta*current(3)
525 absb = sqrt(sum(b(:)**2))
526 bhat(1:ndir) = b(1:ndir) / absb
528 epar = sum(e(:)*bhat(:))
530 call cross(e,bhat,ve)
532 ve(1:ndir) = ve(1:ndir) / absb
533 veabs = sqrt(sum(ve(:)**2))
535 kappa = 1.d0/sqrt(1.0d0 - sum(ve(:)**2)/
c_norm**2)
539 mr = y(ndir+2); upar = y(ndir+1); m=
particle(ipart)%self%m; q=
particle(ipart)%self%q
541 gamma = sqrt(1.0d0+upar**2/
c_norm**2+2.0d0*mr*absb/m/
c_norm**2)*kappa
554 particle(ipart)%payload(1:ndefpayload) = defpayload
568 double precision :: t_s, y(ndir+2)
569 double precision :: dydt(ndir+2)
571 double precision,
dimension(ndir):: vE, b, e, x, bhat, bdotgradb, vEdotgradb, gradkappaB, vfluid, current
572 double precision,
dimension(ndir):: bdotgradvE, vEdotgradvE, u, utmp1, utmp2, utmp3
573 double precision :: upar, Mr, gamma, absb, q, m, epar, kappa
583 x(1:ndir) = y(1:ndir)
588 if (any(x .ne. x))
then
589 write(*,*)
"ERROR IN DERIVS_GCA_RK: NaNs IN X OR Y!"
591 write(*,*)
"y",y(ndir+1:ndir+2)
598 e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) +
particles_eta*current(1)
599 e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) +
particles_eta*current(2)
600 e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) +
particles_eta*current(3)
607 if (any(b .ne. b) .or. any(e .ne. e) &
608 .or. any(bdotgradb .ne. bdotgradb) .or. any(vedotgradb .ne. vedotgradb) &
609 .or. any(gradkappab .ne. gradkappab) .or. any(bdotgradve .ne. bdotgradve) &
610 .or. any(vedotgradve .ne. vedotgradve))
then
611 write(*,*)
"ERROR IN DERIVS_GCA_RK: NaNs IN FIELD QUANTITIES!"
614 write(*,*)
"bdotgradb",bdotgradb
615 write(*,*)
"vEdotgradb",vedotgradb
616 write(*,*)
"gradkappaB",gradkappab
617 write(*,*)
"bdotgradvE",bdotgradve
618 write(*,*)
"vEdotgradvE",vedotgradve
622 absb = sqrt(sum(b(:)**2))
623 if (absb .gt. 0.d0)
then
624 bhat(1:ndir) = b(1:ndir) / absb
628 epar = sum(e(:)*bhat(:))
630 call cross(e,bhat,ve)
631 if (absb .gt. 0.d0) ve(1:ndir) = ve(1:ndir) / absb
634 kappa = 1.d0/sqrt(1.0d0 - sum(ve(:)**2)/
c_norm**2)
635 gamma = sqrt(1.0d0+upar**2/
c_norm**2+2.0d0*mr*absb/m/
c_norm**2)*kappa
641 if (absb .gt. 0.d0)
then
642 utmp1(1:ndir) = bhat(1:ndir)/(absb/kappa**2)
646 utmp2(1:ndir) = mr/(gamma*q)*gradkappab(1:ndir) &
647 + m/q* (upar**2/gamma*bdotgradb(1:ndir) + upar*vedotgradb(1:ndir) &
648 + upar*bdotgradve(1:ndir) + gamma*vedotgradve(1:ndir))
650 utmp2(1:ndir) = utmp2(1:ndir) + upar*epar/(gamma)*ve(1:ndir)
653 call cross(utmp1,utmp2,utmp3)
654 u(1:ndir) = ve(1:ndir) + utmp3(1:ndir)
657 dydt(1:ndir) = ( u(1:ndir) + upar/gamma * bhat(1:ndir) )
658 dydt(ndir+1) = q/m*epar - mr/(m*gamma) * sum(bhat(:)*gradkappab(:)) &
659 + sum(ve(:)*(upar*bdotgradb(:)+gamma*vedotgradb(:)))
667 double precision :: t_s, y(ndir+2)
668 double precision :: dydt(ndir+2)
670 double precision,
dimension(ndir):: vE, b, e, x, bhat, bdotgradb, vEdotgradb, gradkappaB, vfluid, current
671 double precision,
dimension(ndir):: bdotgradvE, vEdotgradvE, u, utmp1, utmp2, utmp3
672 double precision :: upar, Mr, gamma, absb, q, m, epar, kappa
679 x(1:ndir) = y(1:ndir)
684 if (any(x .ne. x))
then
685 call mpistop(
"ERROR IN DERIVS_GCA: NaNs IN X! ABORTING...")
691 e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) +
particles_eta*current(1)
692 e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) +
particles_eta*current(2)
693 e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) +
particles_eta*current(3)
700 absb = sqrt(sum(b(:)**2))
701 if (absb .gt. 0.d0)
then
702 bhat(1:ndir) = b(1:ndir) / absb
707 epar = sum(e(:)*bhat(:))
708 call cross(e,bhat,ve)
709 if (absb .gt. 0.d0) ve(1:ndir) = ve(1:ndir) / absb
712 kappa = sqrt(1.0d0 - sum(ve(:)**2)/
c_norm**2)
713 gamma = sqrt(1.0d0+upar**2/
c_norm**2+2.0d0*mr*absb/m/
c_norm**2)*kappa
718 if (absb .gt. 0.d0)
then
719 utmp1(1:ndir) = bhat(1:ndir)/(absb/kappa**2)
723 utmp2(1:ndir) = mr/(gamma*q)*gradkappab(1:ndir) &
724 + m/q* (upar**2/gamma*bdotgradb(1:ndir) + upar*vedotgradb(1:ndir) &
725 + upar*bdotgradve(1:ndir) + gamma*vedotgradve(1:ndir))
727 utmp2(1:ndir) = utmp2(1:ndir) + upar*epar/(gamma)*ve(1:ndir)
730 call cross(utmp1,utmp2,utmp3)
731 u(1:ndir) = ve(1:ndir) + utmp3(1:ndir)
734 dydt(1:ndir) = ( u(1:ndir) + upar/gamma * bhat(1:ndir) )
735 dydt(ndir+1) = q/m*epar - mr/(m*gamma) * sum(bhat(:)*gradkappab(:)) &
736 + sum(ve(:)*(upar*bdotgradb(:)+gamma*vedotgradb(:)))
742 subroutine gca_update_payload(igrid,w,wold,xgrid,xpart,upart,qpart,mpart,mypayload,mynpayload,particle_time)
744 integer,
intent(in) :: igrid,mynpayload
745 double precision,
intent(in) :: w(ixG^T,1:nw),wold(ixG^T,1:nw)
746 double precision,
intent(in) :: xgrid(ixG^T,1:ndim),xpart(1:ndir),upart(1:ndir),qpart,mpart,particle_time
747 double precision,
intent(out) :: mypayload(mynpayload)
748 double precision,
dimension(1:ndir) :: vE, e, b, bhat, vfluid, current
749 double precision,
dimension(1:ndir) :: drift1, drift2
750 double precision,
dimension(1:ndir) :: drift3, drift4, drift5, drift6, drift7
751 double precision,
dimension(1:ndir) :: bdotgradb, vEdotgradb, gradkappaB
752 double precision,
dimension(1:ndir) :: bdotgradvE, vEdotgradvE
753 double precision,
dimension(1:ndir) :: gradBdrift, reldrift, bdotgradbdrift
754 double precision,
dimension(1:ndir) :: vEdotgradbdrift, bdotgradvEdrift
755 double precision,
dimension(1:ndir) :: vEdotgradvEdrift
756 double precision :: kappa, upar, absb, gamma, vpar, vEabs
757 double precision :: gradBdrift_abs, reldrift_abs, epar
758 double precision :: bdotgradbdrift_abs, vEdotgradbdrift_abs
759 double precision :: bdotgradvEdrift_abs, vEdotgradvEdrift_abs
760 double precision :: momentumpar1, momentumpar2, momentumpar3, momentumpar4
762 call get_vec(
bp, igrid,xpart(1:ndir),particle_time,b)
763 call get_vec(
vp, igrid,xpart(1:ndir),particle_time,vfluid)
764 call get_vec(
jp, igrid,xpart(1:ndir),particle_time,current)
765 e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) +
particles_eta*current(1)
766 e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) +
particles_eta*current(2)
767 e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) +
particles_eta*current(3)
769 absb = sqrt(sum(b(:)**2))
770 bhat(1:ndir) = b(1:ndir) / absb
771 epar = sum(e(:)*bhat(:))
772 call cross(e,bhat,ve)
773 ve(1:ndir) = ve(1:ndir) / absb
774 veabs = sqrt(sum(ve(:)**2))
776 kappa = 1.d0/sqrt(1.0d0 - sum(ve(:)**2)/
c_norm**2)
780 vpar = upart(1)/upart(3)
789 drift1(1:ndir) = bhat(1:ndir)/(absb/kappa**2)
790 drift2(1:ndir) = upart(2)/(upart(3)*qpart)*gradkappab(1:ndir)
792 call cross(drift1,drift2,gradbdrift)
793 gradbdrift_abs = sqrt(sum(gradbdrift(:)**2))
795 drift3(1:ndir) = upar*epar/upart(3)*ve(1:ndir)
796 call cross(drift1,drift3,reldrift)
797 reldrift_abs = sqrt(sum(reldrift(:)**2))
799 drift4(1:ndir) = mpart/qpart* ( upar**2/upart(3)*bdotgradb(1:ndir))
800 call cross(drift1,drift4,bdotgradbdrift)
801 bdotgradbdrift_abs = sqrt(sum(bdotgradbdrift(:)**2))
803 drift5(1:ndir) = mpart/qpart* ( upar*vedotgradb(1:ndir))
804 call cross(drift1,drift5,vedotgradbdrift)
805 vedotgradbdrift_abs = sqrt(sum(vedotgradbdrift(:)**2))
807 drift6(1:ndir) = mpart/qpart* ( upar*bdotgradve(1:ndir))
808 call cross(drift1,drift6,bdotgradvedrift)
809 bdotgradvedrift_abs = sqrt(sum(bdotgradvedrift(:)**2))
811 drift7(1:ndir) = mpart/qpart* (upart(3)*vedotgradve(1:ndir))
812 call cross(drift1,drift7,vedotgradvedrift)
813 vedotgradvedrift_abs = sqrt(sum(vedotgradvedrift(:)**2))
815 momentumpar1 = qpart/mpart*epar
816 momentumpar2 = -(upart(2)/mpart/upart(3))*sum(bhat(:)*gradkappab(:))
817 momentumpar3 = upar*sum(ve(:)*bdotgradb(:))
818 momentumpar4 = upart(3)*sum(ve(:)*vedotgradb(:))
821 if (mynpayload > 0)
then
823 mypayload(1) = sqrt(2.0d0*mpart*upart(2)*absb)/abs(qpart*absb)
825 if (mynpayload > 1)
then
827 mypayload(2) = atan2(sqrt((2.0d0*upart(2)*absb)/(mpart*upart(3)**2)),vpar)
829 if (mynpayload > 2)
then
831 mypayload(3) = sqrt((2.0d0*upart(2)*absb)/(mpart*upart(3)**2))
833 if (mynpayload > 3)
then
835 mypayload(4) = momentumpar1
837 if (mynpayload > 4)
then
839 mypayload(5) = momentumpar2
841 if (mynpayload > 5)
then
843 mypayload(6) = momentumpar3
845 if (mynpayload > 6)
then
847 mypayload(7) = momentumpar4
849 if (mynpayload > 7)
then
853 if (mynpayload > 8)
then
855 mypayload(9) = gradbdrift_abs
857 if (mynpayload > 9)
then
859 mypayload(10) = reldrift_abs
861 if (mynpayload > 10)
then
863 mypayload(11) = bdotgradbdrift_abs
865 if (mynpayload > 11)
then
867 mypayload(12) = vedotgradbdrift_abs
869 if (mynpayload > 12)
then
871 mypayload(13) = bdotgradvedrift_abs
873 if (mynpayload > 13)
then
875 mypayload(14) = vedotgradvedrift_abs
884 double precision,
intent(in) :: end_time
885 double precision :: dt_p
887 double precision :: tout, dt_particles_mype, dt_cfl0, dt_cfl1, dt_a
888 double precision :: dxmin,
vp, a, gammap
890 double precision :: ap0, ap1, dt_cfl_ap0, dt_cfl_ap1, dt_cfl_ap
891 double precision :: dt_euler, dt_tmp
893 double precision :: cfl, uparcfl
894 double precision,
parameter :: uparmin=1.0
d-6*const_c
895 integer :: ipart, iipart, nout, ic^
d, igrid_particle, ipe_particle, ipe
896 logical :: bc_applied
908 dt_tmp = (end_time - partp%self%time)
909 if(dt_tmp .le. 0.0d0) dt_tmp = smalldouble
914 dxmin = min({
rnode(
rpdx^
d_,partp%igrid)},bigdouble)*cfl
917 y(
ndir+1) = partp%self%u(1)
918 y(
ndir+2) = partp%self%u(2)
928 vp = sqrt(sum(v(:)**2))
930 dt_cfl0 = dxmin / max(
vp, smalldouble)
931 dt_cfl_ap0 = uparcfl * abs(max(abs(y(
ndir+1)),uparmin) / max(ap0, smalldouble))
936 dt_euler = min(dt_tmp,dt_cfl0,dt_cfl_ap0)
940 partp%self%u(1) = y(
ndir+1)
941 partp%self%u(2) = y(
ndir+2)
956 vp = sqrt(sum(v(:)**2))
958 dt_cfl1 = dxmin / max(
vp, smalldouble)
959 dt_cfl_ap1 = uparcfl * abs(max(abs(y(
ndir+1)),uparmin) / max(ap1, smalldouble))
962 dt_tmp = min(dt_euler, dt_cfl1, dt_cfl_ap1)
965 partp%self%u(1) = ytmp(
ndir+1)
966 partp%self%u(2) = ytmp(
ndir+2)
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Module with geometry-related routines (e.g., divergence, curl)
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer ixghi
Upper index of grid block arrays.
double precision global_time
The global simulation time.
integer, parameter ndim
Number of spatial dimensions for grid variables.
integer, parameter rpxmin
integer icomm
The MPI communicator.
integer mype
The rank of the current MPI task.
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.
subroutine odeint(ystart, nvar, x1, x2, eps, h1, hmin, nok, nbad, derivs, rkqs, ierror)
subroutine rkqs(y, dydx, n, x, htry, eps, yscal, hdid, hnext, derivs)
Module with shared functionality for all the particle movers.
pure subroutine limit_dt_endtime(t_left, dt_p)
integer num_particles
Number of particles.
double precision const_dt_particles
If positive, a constant time step for the particles.
logical function particle_in_igrid(ipart, igrid)
Quick check if particle is still in igrid.
integer ngridvars
Number of variables for grid field.
integer nusrpayload
Number of user-defined payload variables for a particle.
double precision particles_eta
Resistivity.
subroutine find_particle_ipe(x, igrid_particle, ipe_particle)
pure subroutine get_lfac_from_velocity(v, lfac)
Get Lorentz factor from velocity.
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.
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?
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.
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 gca_update_payload(igrid, w, wold, xgrid, xpart, upart, qpart, mpart, mypayload, mynpayload, particle_time)
Update payload subroutine.
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...
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