MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
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
15contains
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
32 particles_fill_gridvars => advect_fill_gridvars
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')
43 case('ARK4','ARk4','Ark4','ark4')
44 integrator = ark4
45 case default
46 integrator = ark4
47 end select
48
49 particles_integrate => advect_integrate_particles
50
51 end subroutine advect_init
52
54 ! initialise the particles
57
58 double precision :: x(3, num_particles)
59 double precision :: v(3, num_particles)
60 double precision :: q(num_particles)
61 double precision :: m(num_particles)
62 double precision :: rrd(num_particles,ndir)
63 double precision :: w(ixg^t,1:nw)
64 double precision :: defpayload(ndefpayload)
65 double precision :: usrpayload(nusrpayload)
66 integer :: n, idir, igrid, ipe_particle, nparticles_local
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 ! write the first csv file of particles
147
148 end subroutine advect_create_particles
149
150 subroutine advect_fill_gridvars
152
153 double precision, dimension(ixG^T,1:nw) :: w
154 integer :: igrid, iigrid, idir
155
156 ! Fill fluid velocity only
157 do iigrid=1,igridstail; igrid=igrids(iigrid);
158 gridvars(igrid)%w(ixg^t,1:ngridvars) = 0.0d0
159 w(ixg^t,1:nw) = ps(igrid)%w(ixg^t,1:nw)
160 call phys_to_primitive(ixg^ll,ixg^ll,w,ps(igrid)%x)
161 gridvars(igrid)%w(ixg^t,vp(:)) = w(ixg^t,iw_mom(:))
162 gridvars(igrid)%w(ixg^t,rhop) = w(ixg^t,iw_rho)
163 end do
164
165 end subroutine advect_fill_gridvars
166
167 subroutine advect_integrate_particles(end_time)
168 ! this solves dx/dt=v for particles
169 use mod_odeint
172 double precision, intent(in) :: end_time
173
174 double precision, dimension(1:ndir) :: v, x
175 double precision :: defpayload(ndefpayload)
176 double precision :: usrpayload(nusrpayload)
177 double precision :: tloc, tlocnew, dt_p, h1, tk
178 double precision, dimension(1:ndir) :: yk, k1, k2, k3, k4
179 double precision,parameter :: eps=1.0d-6, hmin=1.0d-8
180 integer :: ipart, iipart, igrid
181 integer :: nok, nbad, ierror
182
183 do iipart=1,nparticles_active_on_mype
184 ipart = particles_active_on_mype(iipart);
185 igrid = particle(ipart)%igrid
186 igrid_working = igrid
187 dt_p = advect_get_particle_dt(particle(ipart), end_time)
188 particle(ipart)%self%dt = dt_p
189
190 tloc = particle(ipart)%self%time
191 x(1:ndir) = particle(ipart)%self%x(1:ndir)
192 tlocnew = tloc+dt_p
193
194 ! Position update
195 select case (integrator)
196 case (euler)
197 ! Simple forward Euler
198 call derivs_advect(tloc,x,v)
199 particle(ipart)%self%x(1:ndir) = particle(ipart)%self%x(1:ndir) + dt_p * v(1:ndir)
200
201 case (rk4)
202 ! Runge-Kutta order 4
203 tk = tloc
204 yk = x
205 call derivs_advect(tk,yk,k1)
206 tk = tloc + dt_p/2.d0
207 yk = x + dt_p*k1/2.d0
208 call derivs_advect(tk,yk,k2)
209 tk = tloc + dt_p/2.d0
210 yk = x + dt_p*k2/2.d0
211 call derivs_advect(tk,yk,k3)
212 tk = tloc + dt_p
213 yk = x + dt_p*k3
214 call derivs_advect(tloc,yk,k4)
215 x = x + dt_p/6.d0*(k1 + 2.d0*k2 + 2.d0*k3 + k4)
216
217 case (ark4)
218 ! Adaptive stepwidth RK4:
219 h1 = dt_p/2.0d0
220 call odeint(x,ndir,tloc,tlocnew,eps,h1,hmin,nok,nbad,derivs_advect,rkqs,ierror)
221
222 if (ierror /= 0) then
223 print *, "odeint returned error code", ierror
224 print *, "1 means hmin too small, 2 means MAXSTP exceeded"
225 print *, "Having a problem with particle", iipart
226 end if
227 end select
228 particle(ipart)%self%x(1:ndir) = x(1:ndir)
229
230 ! Velocity update
231 call get_vec_advect(igrid,x,tlocnew,v,vp(1),vp(ndir))
232 particle(ipart)%self%u(1:ndir) = v(1:ndir)
233
234 ! Time update
235 particle(ipart)%self%time = tlocnew
236
237 ! Update payload
238 if(ndefpayload>0) then
239 call advect_update_payload(igrid,x,v,0.d0,0.d0,defpayload,ndefpayload,tlocnew)
240 particle(ipart)%payload(1:ndefpayload) = defpayload
241 end if
242 if (associated(usr_update_payload)) then
243 call usr_update_payload(igrid,x,v,0.d0,0.d0,usrpayload,nusrpayload,tlocnew)
244 particle(ipart)%payload(ndefpayload+1:npayload) = usrpayload
245 end if
246
247 end do
248
249 end subroutine advect_integrate_particles
250
251 !> Payload update
252 subroutine advect_update_payload(igrid,xpart,upart,qpart,mpart,mypayload,mynpayload,particle_time)
254 integer, intent(in) :: igrid,mynpayload
255 double precision, intent(in) :: xpart(1:ndir),upart(1:ndir),qpart,mpart,particle_time
256 double precision, intent(out) :: mypayload(mynpayload)
257 double precision :: rho, rho1, rho2, td
258
259 td = (particle_time - global_time) / dt
260
261 ! Payload 1 is density
262 if (mynpayload > 0 ) then
263 if (time_advance) then
264 call interpolate_var(igrid,ixg^ll,ixm^ll,gridvars(igrid)%wold(ixg^t,rhop),ps(igrid)%x,xpart,rho1)
265 call interpolate_var(igrid,ixg^ll,ixm^ll,gridvars(igrid)%w(ixg^t,rhop),ps(igrid)%x,xpart,rho2)
266 rho = rho1 * (1.0d0 - td) + rho2 * td
267 else
268 call interpolate_var(igrid,ixg^ll,ixm^ll,gridvars(igrid)%w(ixg^t,rhop),ps(igrid)%x,xpart,rho)
269 end if
270 mypayload(1) = rho * w_convert_factor(1)
271 end if
272
273 end subroutine advect_update_payload
274
275 subroutine derivs_advect(t_s,x,dxdt)
277 use mod_geometry
278 double precision :: t_s, x(ndir)
279 double precision :: dxdt(ndir)
280
281 double precision :: v(ndir)
282
283 call get_vec_advect(igrid_working,x,t_s,v,vp(1),vp(ndir))
284 select case (coordinate)
285 case (cylindrical)
286 v(phi_) = v(phi_)/x(r_)
287 case (spherical)
288 v(2) = v(2)/x(1)
289 v(3) = v(3)/(x(1)*sin(x(2)))
290 end select
291 dxdt(:) = v(:)
292
293 end subroutine derivs_advect
294
295 function advect_get_particle_dt(partp, end_time) result(dt_p)
297 type(particle_ptr), intent(in) :: partp
298 double precision, intent(in) :: end_time
299 double precision :: dt_p
300 double precision :: v(1:ndir),dtdims(1:ndim)
301 integer :: ipart, iipart, idims
302
303 if (const_dt_particles > 0) then
304 dt_p = const_dt_particles
305 return
306 end if
307
308 ! make sure we step only one cell at a time:
309 call derivs_advect(partp%self%time,partp%self%x,v)
310 do idims=1,ndim
311 dtdims(idims)=minval(ps(partp%igrid)%ds(ixm^t,idims))/(max(abs(v(idims)),smalldouble))
312 end do
313
314 dt_p = particles_cfl*minval(dtdims)
315
316 ! Make sure we do not advance beyond end_time
317 call limit_dt_endtime(end_time - partp%self%time, dt_p)
318
319 end function advect_get_particle_dt
320
321 subroutine get_vec_advect(igrid,x,tloc,var,ibeg,iend)
323
324 integer,intent(in) :: igrid, ibeg, iend
325 double precision,dimension(ndir), intent(in) :: x
326 double precision, intent(in) :: tloc
327 double precision,dimension(iend-ibeg+1), intent(out) :: var
328 double precision,dimension(iend-ibeg+1) :: e1, e2
329 double precision :: td
330 integer :: ivar, iloc
331
332 if(.not.time_advance) then
333 do ivar=ibeg,iend
334 iloc = ivar-ibeg+1
335 call interpolate_var(igrid,ixg^ll,ixm^ll,gridvars(igrid)%w(ixg^t,ivar),ps(igrid)%x(ixg^t,1:ndim),x,var(iloc))
336 end do
337 else
338 td = (tloc - global_time) / dt
339 do ivar=ibeg,iend
340 iloc = ivar-ibeg+1
341 call interpolate_var(igrid,ixg^ll,ixm^ll,gridvars(igrid)%wold(ixg^t,ivar),ps(igrid)%x(ixg^t,1:ndim),x,e1(iloc))
342 call interpolate_var(igrid,ixg^ll,ixm^ll,gridvars(igrid)%w(ixg^t,ivar),ps(igrid)%x(ixg^t,1:ndim),x,e2(iloc))
343 var(iloc) = e1(iloc) * (1.0d0 - td) + e2(iloc) * td
344 end do
345 end if
346
347 end subroutine get_vec_advect
348
349end 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
integer, parameter cylindrical
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.
double precision, dimension(:), allocatable, parameter d
double precision dt
global time step
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.
logical time_advance
do time evolving
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, 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.
subroutine write_particle_output()
double precision dtsave_particles
Time interval to save 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
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:
subroutine interpolate_var(igrid, ixil, ixol, gf, x, xloc, gfloc)
double precision particles_cfl
Particle CFL safety factor.
integer rhop
Variable index for density.
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