Particles II – bhac.science

Custom initialization, injection and destruction or particles

The routines that control the initialization, injection, and destruction of particles are located in the file mod_particles_user.t in $BHAC_DIR/src/particles/ . To customize this functions, we copy this file to our working directory.

cp $BHAC_DIR/src/particles/mod_particles_user.t ./

As with all of the source code, the copy of the file in our working directory will override that in $BHAC_DIR.

Particles are initialized in the routine init_particles_user() . The default version of this routine initializes particles at random positions until the pre-defined number of particles is reached. This is done in the loop from lines 82 to 107.

82 do while (nparticles_added .lt. nparticles_init)
83 nparticles_added = nparticles_added + 1
84
85 {^D&x(^D) = xprobmin^D + r^D(nparticles_added) * (xprobmax^D - xprobmin^D)}
…

If the particle is in our CPU, its attributes are assigned in lines 92-103

92 particles_to_inject(ninject)%self%follow = follow(nparticles_added)
93 particles_to_inject(ninject)%self%q = - CONST_e
94 particles_to_inject(ninject)%self%m = CONST_me
95
96 particles_to_inject(ninject)%self%t = t
97 particles_to_inject(ninject)%self%dt = 0.0d0
98
99 particles_to_inject(ninject)%self%x(:) = x
100
101 call get_vec(igrid_particle,x,0.0d0,v,vp1_,vp^NC_)
102 particles_to_inject(ninject)%self%u(:) = v(:)
103 particles_to_inject(ninject)%payload(:) = 666
104

Here we will edit the function so that it initializes particles based on the local density. The FMtorus set-up normalizes the maximum density to 1, so we can safely use the density ‘rho’ as the probability to initialize a particle.

We need to add a few new variables:

! … local …
double precision :: local_rho,rand,r^D

and change slightly the structure of the loop so that every time a position is generated randomly, it also generates a random number rand between 0 and 1 and compares it with the local density, local_rho, placing a particle if rand<local_rho.

For this, we need to interpolate the local density to the particle location, which we can accomplish using the function function interpolate_var(). The content of the full routine now looks:

!=============================================================================
subroutine initial_particles_user()
! add the initial particles
use constants
use mod_random
include 'amrvacdef.f90'

double precision, dimension(ndir)       :: x
integer                                 :: igrid_particle, ipe_particle
integer(i8)                             :: seed(2)
type(rng_t)                             :: myrand
double precision                        :: v(1:ndir)
logical, dimension(:),allocatable       :: follow
type(particle_ptr), dimension(:), allocatable    :: particles_to_inject
integer                                 :: ninject, nparticles_added, ipart
! ... local ...
double precision :: local_rho,rand,r^D
!-----------------------------------------------------------------------------

allocate(follow(nparticles_init), particles_to_inject(nparticles_init))
do ipart=1,nparticles_init
   allocate(particles_to_inject(ipart)%self)
   allocate(particles_to_inject(ipart)%payload(npayload))
end do

! initialise the random number generator
seed = [310952_i8,24378_i8]
call myrand%set_seed(seed)

! flags to follow particles
follow      = .false.
follow(1:3)  = .true.


ninject = 0 
nparticles_added = 0 
x=0 
do while (nparticles_added .lt. nparticles_init)
   {^D& r^D=myrand%unif_01()\}

   {^D&x(^D) = xprobmin^D + r^D * (xprobmax^D - xprobmin^D)\}

   call find_particle_ipe(x,igrid_particle,ipe_particle)
   ! Interpolate local density from grid
   call interpolate_var(igrid_particle,ixG^LL,ixM^LL,pw(igrid_particle)%w(ixG^T,rho_),&
                        px(igrid_particle)%x(ixG^T,1:ndim),x,local_rho)

   rand=myrand%unif_01()
   print *, rand,local_rho
   if (rand.lt.local_rho) then
       print *, 'Particle added'
       nparticles_added = nparticles_added + 1 

       if (ipe_particle == mype) then 
          ninject = ninject + 1 

          particles_to_inject(ninject)%self%follow = follow(nparticles_added)
          particles_to_inject(ninject)%self%q      = - CONST_e
          particles_to_inject(ninject)%self%m      =   CONST_me

          particles_to_inject(ninject)%self%t      = t 
          particles_to_inject(ninject)%self%dt     = 0.0d0

          particles_to_inject(ninject)%self%x(:)   = x 

          call get_vec(igrid_particle,x,0.0d0,v,vp1_,vp^NC_)
          particles_to_inject(ninject)%self%u(:) = v(:)
          particles_to_inject(ninject)%payload(:) = 666 

       end if
   end if

end do

call inject_particles(particles_to_inject(1:ninject))
end subroutine initial_particles_user
!=============================================================================

The following plot shows the location of particles at the initial snapshot before and after modifying this routine. The routines for particle injection add_particles_user and particle destruction destroy_user can be modified in a similar manner.

Custom initialization, injection and destruction or particles

Payloads