MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_particle_advect.t
Go to the documentation of this file.
1 !> Tracer for advected particles moving with fluid flows
2 !> By Jannis Teunissen, Bart Ripperda, Oliver Porth, and Fabio Bacchini (2017-2020)
5 
6  private
7 
8  public :: advect_init
10  integer, parameter :: Euler = 1, rk4=2, ark4=3
11 
12  ! Variables
13  public :: vp, rhop
14 
15 contains
16 
17  subroutine advect_init()
19  integer :: idir
20 
21  ngridvars=0
22  allocate(vp(ndir))
23  do idir = 1, ndir
24  vp(idir) = idir
25  end do
27  if(ndefpayload>0) then
28  rhop=vp(ndir)+1
30  end if
31 
33 
34  if (associated(particles_define_additional_gridvars)) then
36  end if
37 
38  select case(integrator_type_particles)
39  case('Euler','euler')
40  integrator = euler
41  case('RK4','Rk4','rk4')
42  integrator = rk4
43  case('ARK4','ARk4','Ark4','ark4')
44  integrator = ark4
45  case default
46  integrator = ark4
47  end select
48 
50 
51  end subroutine advect_init
52 
54  ! initialise the particles
57 
58  integer :: n, idir, igrid, ipe_particle, nparticles_local
59  double precision :: x(3, num_particles)
60  double precision :: v(3, num_particles)
61  double precision :: q(num_particles)
62  double precision :: m(num_particles)
63  double precision :: rrd(num_particles,ndir)
64  double precision :: w(ixg^t,1:nw)
65  double precision :: defpayload(ndefpayload)
66  double precision :: usrpayload(nusrpayload)
67  logical :: follow(num_particles), check
68 
69  follow = .false.
70  x = 0.0d0
71 
72  if (mype==0) then
73  if (.not. associated(usr_create_particles)) then
74  ! Randomly distributed
75  do idir=1,ndir
76  do n = 1, num_particles
77  rrd(n,idir) = rng%unif_01()
78  end do
79  end do
80  do n=1, num_particles
81  {^d&x(^d,n) = xprobmin^d + rrd(n,^d) * (xprobmax^d - xprobmin^d)\}
82  end do
83  else
84  call usr_create_particles(num_particles, x, v, q, m, follow)
85  end if
86  end if
87 
88  call mpi_bcast(x,3*num_particles,mpi_double_precision,0,icomm,ierrmpi)
89  call mpi_bcast(follow,num_particles,mpi_logical,0,icomm,ierrmpi)
90 
91  nparticles_local = 0
92 
93  do n=1,num_particles
94  call find_particle_ipe(x(:,n),igrid,ipe_particle)
95  particle(n)%igrid = igrid
96  particle(n)%ipe = ipe_particle
97 
98  if(ipe_particle == mype) then
99  check = .true.
100 
101  ! Check for user-defined modifications or rejection conditions
102  if (associated(usr_check_particle)) call usr_check_particle(igrid, x(:,n), v(:,n), q(n), m(n), follow(n), check)
103  if (check) then
105  else
106  cycle
107  end if
108 
109  nparticles_local = nparticles_local + 1
110 
111  allocate(particle(n)%self)
112  particle(n)%self%follow = follow(n)
113  particle(n)%self%index = n
114  particle(n)%self%time = global_time
115  particle(n)%self%dt = 0.0d0
116  particle(n)%self%x = 0.d0
117  particle(n)%self%x(:) = x(:,n)
118  w=ps(igrid)%w
119  call phys_to_primitive(ixg^ll,ixg^ll,w,ps(igrid)%x)
120  do idir=1,ndir
121  call interpolate_var(igrid,ixg^ll,ixm^ll,&
122  w(ixg^t,iw_mom(idir)),ps(igrid)%x,x(:,n),v(idir,n))
123  end do
124  particle(n)%self%u(:) = 0.d0
125  particle(n)%self%u(1:ndir) = v(1:ndir,n)
126 
127  ! Compute default and user-defined payloads
128  allocate(particle(n)%payload(npayload))
129  if(ndefpayload>0) then
130  call advect_update_payload(igrid,x(:,n),v(:,n),q(n),m(n),defpayload,ndefpayload,0.d0)
131  particle(n)%payload(1:ndefpayload)=defpayload
132  end if
133  if (associated(usr_update_payload)) then
134  call usr_update_payload(igrid,x(:,n),v(:,n),q(n),m(n),usrpayload,nusrpayload,0.d0)
135  particle(n)%payload(ndefpayload+1:npayload)=usrpayload
136  end if
137  end if
138 
139  end do
140 
141  call mpi_allreduce(nparticles_local,nparticles,1,mpi_integer,mpi_sum,icomm,ierrmpi)
142 
143  end subroutine advect_create_particles
144 
147 
148  integer :: igrid, iigrid, idir
149  double precision, dimension(ixG^T,1:nw) :: w
150 
151  ! Fill fluid velocity only
152  do iigrid=1,igridstail; igrid=igrids(iigrid);
153  gridvars(igrid)%w(ixg^t,1:ngridvars) = 0.0d0
154  w(ixg^t,1:nw) = ps(igrid)%w(ixg^t,1:nw)
155  call phys_to_primitive(ixg^ll,ixg^ll,w,ps(igrid)%x)
156  gridvars(igrid)%w(ixg^t,vp(:)) = w(ixg^t,iw_mom(:))
157  gridvars(igrid)%w(ixg^t,rhop) = w(ixg^t,iw_rho)
158  end do
159 
160  end subroutine advect_fill_gridvars
161 
162  subroutine advect_integrate_particles(end_time)
163  ! this solves dx/dt=v for particles
164  use mod_odeint
167  double precision, intent(in) :: end_time
168 
169  double precision, dimension(1:ndir) :: v, x
170  double precision :: defpayload(ndefpayload)
171  double precision :: usrpayload(nusrpayload)
172  double precision :: tloc, tlocnew, dt_p, h1, tk
173  double precision, dimension(1:ndir) :: yk, k1, k2, k3, k4
174  double precision,parameter :: eps=1.0d-6, hmin=1.0d-8
175  integer :: ipart, iipart, igrid
176  integer :: nok, nbad, ierror
177 
178  do iipart=1,nparticles_active_on_mype
179  ipart = particles_active_on_mype(iipart);
180  igrid = particle(ipart)%igrid
181  igrid_working = igrid
182  dt_p = advect_get_particle_dt(particle(ipart), end_time)
183  particle(ipart)%self%dt = dt_p
184 
185  tloc = particle(ipart)%self%time
186  x(1:ndir) = particle(ipart)%self%x(1:ndir)
187  tlocnew = tloc+dt_p
188 
189  ! Position update
190  select case (integrator)
191  case (euler)
192  ! Simple forward Euler
193  call derivs_advect(tloc,x,v)
194  particle(ipart)%self%x(1:ndir) = particle(ipart)%self%x(1:ndir) + dt_p * v(1:ndir)
195 
196  case (rk4)
197  ! Runge-Kutta order 4
198  tk = tloc
199  yk = x
200  call derivs_advect(tk,yk,k1)
201  tk = tloc + dt_p/2.d0
202  yk = x + dt_p*k1/2.d0
203  call derivs_advect(tk,yk,k2)
204  tk = tloc + dt_p/2.d0
205  yk = x + dt_p*k2/2.d0
206  call derivs_advect(tk,yk,k3)
207  tk = tloc + dt_p
208  yk = x + dt_p*k3
209  call derivs_advect(tloc,yk,k4)
210  x = x + dt_p/6.d0*(k1 + 2.d0*k2 + 2.d0*k3 + k4)
211 
212  case (ark4)
213  ! Adaptive stepwidth RK4:
214  h1 = dt_p/2.0d0
215  call odeint(x,ndir,tloc,tlocnew,eps,h1,hmin,nok,nbad,derivs_advect,rkqs,ierror)
216 
217  if (ierror /= 0) then
218  print *, "odeint returned error code", ierror
219  print *, "1 means hmin too small, 2 means MAXSTP exceeded"
220  print *, "Having a problem with particle", iipart
221  end if
222  end select
223  particle(ipart)%self%x(1:ndir) = x(1:ndir)
224 
225  ! Velocity update
226  call get_vec_advect(igrid,x,tlocnew,v,vp(1),vp(ndir))
227  particle(ipart)%self%u(1:ndir) = v(1:ndir)
228 
229  ! Time update
230  particle(ipart)%self%time = tlocnew
231 
232  ! Update payload
233  if(ndefpayload>0) then
234  call advect_update_payload(igrid,x,v,0.d0,0.d0,defpayload,ndefpayload,tlocnew)
235  particle(ipart)%payload(1:ndefpayload) = defpayload
236  end if
237  if (associated(usr_update_payload)) then
238  call usr_update_payload(igrid,x,v,0.d0,0.d0,usrpayload,nusrpayload,tlocnew)
239  particle(ipart)%payload(ndefpayload+1:npayload) = usrpayload
240  end if
241 
242  end do
243 
244  end subroutine advect_integrate_particles
245 
246  !> Payload update
247  subroutine advect_update_payload(igrid,xpart,upart,qpart,mpart,mypayload,mynpayload,particle_time)
249  integer, intent(in) :: igrid,mynpayload
250  double precision, intent(in) :: xpart(1:ndir),upart(1:ndir),qpart,mpart,particle_time
251  double precision, intent(out) :: mypayload(mynpayload)
252  double precision :: rho, rho1, rho2, td
253 
254  td = (particle_time - global_time) / dt
255 
256  ! Payload 1 is density
257  if (mynpayload > 0 ) then
258  if (time_advance) then
259  call interpolate_var(igrid,ixg^ll,ixm^ll,gridvars(igrid)%wold(ixg^t,rhop),ps(igrid)%x,xpart,rho1)
260  call interpolate_var(igrid,ixg^ll,ixm^ll,gridvars(igrid)%w(ixg^t,rhop),ps(igrid)%x,xpart,rho2)
261  rho = rho1 * (1.0d0 - td) + rho2 * td
262  else
263  call interpolate_var(igrid,ixg^ll,ixm^ll,gridvars(igrid)%w(ixg^t,rhop),ps(igrid)%x,xpart,rho)
264  end if
265  mypayload(1) = rho * w_convert_factor(1)
266  end if
267 
268  end subroutine advect_update_payload
269 
270  subroutine derivs_advect(t_s,x,dxdt)
272  use mod_geometry
273  double precision :: t_s, x(ndir)
274  double precision :: dxdt(ndir)
275 
276  double precision :: v(ndir)
277 
278  call get_vec_advect(igrid_working,x,t_s,v,vp(1),vp(ndir))
279  select case (coordinate)
280  case (cylindrical)
281  v(phi_) = v(phi_)/x(r_)
282  case (spherical)
283  v(2) = v(2)/x(1)
284  v(3) = v(3)/(x(1)*sin(x(2)))
285  end select
286  dxdt(:) = v(:)
287 
288  end subroutine derivs_advect
289 
290  function advect_get_particle_dt(partp, end_time) result(dt_p)
292  type(particle_ptr), intent(in) :: partp
293  double precision, intent(in) :: end_time
294  double precision :: dt_p
295  integer :: ipart, iipart, idims
296  double precision :: v(1:ndir),dtdims(1:ndim)
297 
298  if (const_dt_particles > 0) then
299  dt_p = const_dt_particles
300  return
301  end if
302 
303  ! make sure we step only one cell at a time:
304  call derivs_advect(partp%self%time,partp%self%x,v)
305  do idims=1,ndim
306  dtdims(idims)=minval(ps(partp%igrid)%ds(ixm^t,idims))/(max(abs(v(idims)),smalldouble))
307  end do
308 
309  dt_p = particles_cfl*minval(dtdims)
310 
311  ! Make sure we do not advance beyond end_time
312  call limit_dt_endtime(end_time - partp%self%time, dt_p)
313 
314  end function advect_get_particle_dt
315 
316  subroutine get_vec_advect(igrid,x,tloc,var,ibeg,iend)
318 
319  integer,intent(in) :: igrid, ibeg, iend
320  double precision,dimension(ndir), intent(in) :: x
321  double precision, intent(in) :: tloc
322  double precision,dimension(iend-ibeg+1), intent(out) :: var
323  double precision,dimension(iend-ibeg+1) :: e1, e2
324  integer :: ivar, iloc
325  double precision :: td
326 
327  if(.not.time_advance) then
328  do ivar=ibeg,iend
329  iloc = ivar-ibeg+1
330  call interpolate_var(igrid,ixg^ll,ixm^ll,gridvars(igrid)%w(ixg^t,ivar),ps(igrid)%x(ixg^t,1:ndim),x,var(iloc))
331  end do
332  else
333  td = (tloc - global_time) / dt
334  do ivar=ibeg,iend
335  iloc = ivar-ibeg+1
336  call interpolate_var(igrid,ixg^ll,ixm^ll,gridvars(igrid)%wold(ixg^t,ivar),ps(igrid)%x(ixg^t,1:ndim),x,e1(iloc))
337  call interpolate_var(igrid,ixg^ll,ixm^ll,gridvars(igrid)%w(ixg^t,ivar),ps(igrid)%x(ixg^t,1:ndim),x,e2(iloc))
338  var(iloc) = e1(iloc) * (1.0d0 - td) + e2(iloc) * td
339  end do
340  end if
341 
342  end subroutine get_vec_advect
343 
344 end module mod_particle_advect
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
integer coordinate
Definition: mod_geometry.t:7
integer, parameter spherical
Definition: mod_geometry.t:11
integer, parameter cylindrical
Definition: mod_geometry.t:10
This module contains definitions of global parameters and variables and some generic functions/subrou...
double precision, dimension(:), allocatable w_convert_factor
Conversion factors the primitive variables.
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 ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
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
Tracer for advected particles moving with fluid flows By Jannis Teunissen, Bart Ripperda,...
subroutine, public advect_init()
subroutine advect_integrate_particles(end_time)
subroutine derivs_advect(t_s, x, dxdt)
subroutine advect_update_payload(igrid, xpart, upart, qpart, mpart, mypayload, mynpayload, particle_time)
Payload update.
double precision function advect_get_particle_dt(partp, end_time)
subroutine get_vec_advect(igrid, x, tloc, var, ibeg, iend)
subroutine, public advect_create_particles()
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.
integer ngridvars
Number of variables for grid field.
integer nusrpayload
Number of user-defined payload variables for a particle.
subroutine find_particle_ipe(x, igrid_particle, ipe_particle)
integer integrator
Integrator to be used for particles.
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 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.
type(particle_ptr), dimension(:), allocatable particle
integer, dimension(:), allocatable vp
Variable index for fluid velocity.
procedure(sub_integrate), pointer particles_integrate
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:
double precision particles_cfl
Particle CFL safety factor.
integer rhop
Variable index for density.
subroutine interpolate_var(igrid, ixIL, ixOL, gf, x, xloc, gfloc)
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