10 integer,
parameter :: Euler = 1, rk4=2, ark4=3
41 case(
'RK4',
'Rk4',
'rk4')
43 case(
'ARK4',
'ARk4',
'Ark4',
'ark4')
63 double precision :: w(ixg^t,1:nw)
66 integer :: n, idir, igrid, ipe_particle, nparticles_local
77 rrd(n,idir) =
rng%unif_01()
81 {^
d&x(^
d,n) = xprobmin^
d + rrd(n,^
d) * (xprobmax^
d - xprobmin^
d)\}
98 if(ipe_particle ==
mype)
then
109 nparticles_local = nparticles_local + 1
119 call phys_to_primitive(ixg^
ll,ixg^
ll,w,ps(igrid)%x)
122 w(ixg^t,iw_mom(idir)),ps(igrid)%x,x(:,n),v(idir,n))
130 call advect_update_payload(igrid,x(:,n),v(:,n),q(n),m(n),defpayload,
ndefpayload,0.d0)
150 subroutine advect_fill_gridvars
153 double precision,
dimension(ixG^T,1:nw) :: w
154 integer :: igrid, iigrid, idir
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)
165 end subroutine advect_fill_gridvars
167 subroutine advect_integrate_particles(end_time)
172 double precision,
intent(in) :: end_time
174 double precision,
dimension(1:ndir) :: v, x
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.0
d-6, hmin=1.0
d-8
180 integer :: ipart, iipart, igrid
181 integer :: nok, nbad, ierror
187 dt_p = advect_get_particle_dt(
particle(ipart), end_time)
198 call derivs_advect(tloc,x,v)
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)
214 call derivs_advect(tloc,yk,k4)
215 x = x + dt_p/6.d0*(k1 + 2.d0*k2 + 2.d0*k3 + k4)
220 call odeint(x,
ndir,tloc,tlocnew,eps,h1,hmin,nok,nbad,derivs_advect,
rkqs,ierror)
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
231 call get_vec_advect(igrid,x,tlocnew,v,
vp(1),
vp(
ndir))
239 call advect_update_payload(igrid,x,v,0.d0,0.d0,defpayload,
ndefpayload,tlocnew)
249 end subroutine advect_integrate_particles
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
262 if (mynpayload > 0 )
then
266 rho = rho1 * (1.0d0 - td) + rho2 * td
273 end subroutine advect_update_payload
275 subroutine derivs_advect(t_s,x,dxdt)
278 double precision :: t_s, x(
ndir)
279 double precision :: dxdt(
ndir)
281 double precision :: v(
ndir)
289 v(3) = v(3)/(x(1)*sin(x(2)))
293 end subroutine derivs_advect
295 function advect_get_particle_dt(partp, end_time)
result(dt_p)
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
309 call derivs_advect(partp%self%time,partp%self%x,v)
311 dtdims(idims)=minval(ps(partp%igrid)%ds(
ixm^t,idims))/(max(abs(v(idims)),smalldouble))
319 end function advect_get_particle_dt
321 subroutine get_vec_advect(igrid,x,tloc,var,ibeg,iend)
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
343 var(iloc) = e1(iloc) * (1.0d0 - td) + e2(iloc) * td
347 end subroutine get_vec_advect
Module with geometry-related routines (e.g., divergence, curl)
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.
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)
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