MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_particle_sample.t
Go to the documentation of this file.
1 !> Scattered sampling based on fixed- or moving-particle interpolation
2 !> By Fabio Bacchini (2020)
5 
6  private
7 
8  public :: sample_init
10 
11 contains
12 
13  subroutine sample_init()
15  integer :: idir
16 
17  ngridvars=nw
18 
20 
21  if (associated(particles_define_additional_gridvars)) then
23  end if
24 
26 
27  end subroutine sample_init
28 
30  ! initialise the particles (=fixed interpolation points)
34 
35  integer :: n, idir, igrid, ipe_particle, nparticles_local
36  double precision :: x(3, num_particles)
37  double precision :: v(3, num_particles)
38  double precision :: q(num_particles)
39  double precision :: m(num_particles)
40  double precision :: rrd(num_particles,ndir)
41  double precision :: w(ixg^t,1:nw)
42  double precision :: defpayload(ndefpayload)
43  double precision :: usrpayload(nusrpayload)
44  logical :: follow(num_particles), check
45 
46  follow = .false.
47  x = 0.0d0
48 
49  if (mype==0) then
50  if (.not. associated(usr_create_particles)) then
51  ! Randomly distributed
52  do idir=1,ndir
53  do n = 1, num_particles
54  rrd(n,idir) = rng%unif_01()
55  end do
56  end do
57  do n=1, num_particles
58  {^d&x(^d,n) = xprobmin^d + rrd(n+1,^d) * (xprobmax^d - xprobmin^d)\}
59  end do
60  else
61  call usr_create_particles(num_particles, x, v, q, m, follow)
62  end if
63  end if
64 
65  call mpi_bcast(x,3*num_particles,mpi_double_precision,0,icomm,ierrmpi)
66  call mpi_bcast(follow,num_particles,mpi_logical,0,icomm,ierrmpi)
67 
68  nparticles_local = 0
69 
70  do n=1,num_particles
71  call find_particle_ipe(x(:,n),igrid,ipe_particle)
72  particle(n)%igrid = igrid
73  particle(n)%ipe = ipe_particle
74 
75  if(ipe_particle == mype) then
76  check = .true.
77 
78  ! Check for user-defined modifications or rejection conditions
79  if (associated(usr_check_particle)) call usr_check_particle(igrid, x(:,n), v(:,n), q(n), m(n), follow(n), check)
80  if (check) then
82  else
83  cycle
84  end if
85 
86  nparticles_local = nparticles_local + 1
87 
88  allocate(particle(n)%self)
89  particle(n)%self%follow = follow(n)
90  particle(n)%self%index = n
91  particle(n)%self%time = global_time
92  particle(n)%self%dt = 0.0d0
93  particle(n)%self%x = 0.d0
94  particle(n)%self%x(:) = x(:,n)
95  particle(n)%self%u(:) = 0.d0
96 
97  allocate(particle(n)%payload(npayload))
98  call sample_update_payload(igrid,x(:,n),v(:,n),q(n),m(n),defpayload,ndefpayload,0.d0)
99  particle(n)%payload(1:ndefpayload) = defpayload
100  if (associated(usr_update_payload)) then
101  call usr_update_payload(igrid,x(:,n),v(:,n),q(n),m(n),usrpayload,nusrpayload,0.d0)
102  particle(n)%payload(ndefpayload+1:npayload)=usrpayload
103  end if
104  end if
105 
106  end do
107 
108  call mpi_allreduce(nparticles_local,nparticles,1,mpi_integer,mpi_sum,icomm,ierrmpi)
109 
110  end subroutine sample_create_particles
111 
114 
115  integer :: igrid, iigrid, idir
116  double precision, dimension(ixG^T,1:nw) :: w
117 
118  do iigrid=1,igridstail; igrid=igrids(iigrid);
119 
120  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
121  w(ixg^t,1:nw) = ps(igrid)%w(ixg^t,1:nw)
122  call phys_to_primitive(ixg^ll,ixg^ll,w,ps(igrid)%x)
123  ! fill all variables:
124  gridvars(igrid)%w(ixg^t,1:ngridvars) = w(ixg^t,1:ngridvars)
125 
126  end do
127 
128  end subroutine sample_fill_gridvars
129 
130  subroutine sample_integrate_particles(end_time)
131  ! this interpolates the HD/MHD quantities at the particle positions
134  double precision, intent(in) :: end_time
135 
136  double precision, dimension(1:ndir) :: v, x
137  double precision :: defpayload(ndefpayload)
138  double precision :: usrpayload(nusrpayload)
139  double precision :: tloc, tlocnew, dt_p, h1
140  double precision,parameter :: eps=1.0d-6, hmin=1.0d-8
141  integer :: ipart, iipart, igrid
142  integer :: nok, nbad, ierror
143 
144  do iipart=1,nparticles_active_on_mype
145  ipart = particles_active_on_mype(iipart);
146  dt_p = sample_get_particle_dt(particle(ipart), end_time)
147  particle(ipart)%self%dt = dt_p
148 
149  igrid = particle(ipart)%igrid
150  igrid_working = igrid
151  tloc = particle(ipart)%self%time
152  x(1:ndir) = particle(ipart)%self%x(1:ndir)
153  tlocnew = tloc+dt_p
154 
155  ! Position update (if defined)
156  ! TODO: this may create problems with interpolation out of boundaries
157  if (associated(usr_particle_position)) call usr_particle_position(x,particle(ipart)%self%index,tloc,tlocnew)
158  particle(ipart)%self%x(1:ndir) = x
159 
160  ! Time update
161  particle(ipart)%self%time = tlocnew
162 
163  ! Update payload
164  call sample_update_payload(igrid,x,v,0.d0,0.d0,defpayload,ndefpayload,tlocnew)
165  particle(ipart)%payload(1:ndefpayload) = defpayload
166  if (associated(usr_update_payload)) then
167  call usr_update_payload(igrid,x,v,0.d0,0.d0,usrpayload,nusrpayload,tlocnew)
168  particle(ipart)%payload(ndefpayload+1:npayload) = usrpayload
169  end if
170 
171  end do
172 
173  end subroutine sample_integrate_particles
174 
175  !> Payload update
176  subroutine sample_update_payload(igrid,xpart,upart,qpart,mpart,mypayload,mynpayload,particle_time)
178  integer, intent(in) :: igrid,mynpayload
179  double precision, intent(in) :: xpart(1:ndir),upart(1:ndir),qpart,mpart,particle_time
180  double precision, intent(out) :: mypayload(mynpayload)
181  double precision :: myw(ixG^T,1:nw),mywold(ixG^T,1:nw)
182  double precision :: wp, wpold, td
183  integer :: ii
184 
185 
186  ! There are npayload=nw payloads, one for each primitive fluid quantity
187  myw(ixg^t,1:nw) = gridvars(igrid)%w(ixg^t,1:nw)
188  if (time_advance) mywold(ixg^t,1:nw) = gridvars(igrid)%wold(ixg^t,1:nw)
189 
190  if (.not.saveprim) then
191  call phys_to_conserved(ixg^ll,ixg^ll,myw,ps(igrid)%x)
192  if (time_advance) call phys_to_conserved(ixg^ll,ixg^ll,mywold,ps(igrid)%x)
193  end if
194 
195  do ii=1,mynpayload
196  call interpolate_var(igrid,ixg^ll,ixm^ll,myw(ixg^t,ii),ps(igrid)%x,xpart,wp)
197  if (time_advance) then
198  td = (particle_time - global_time) / dt
199  call interpolate_var(igrid,ixg^ll,ixm^ll,mywold(ixg^t,ii),ps(igrid)%x,xpart,wpold)
200  wp = wpold * (1.0d0 - td) + wp * td
201  end if
202  mypayload(ii) = wp*w_convert_factor(ii)
203  end do
204 
205  end subroutine sample_update_payload
206 
207  function sample_get_particle_dt(partp, end_time) result(dt_p)
209  use mod_geometry
211  type(particle_ptr), intent(in) :: partp
212  double precision, intent(in) :: end_time
213  double precision :: dt_p
214  integer :: ipart, iipart, nout, id
215  double precision :: tout, dt_cfl
216  double precision :: v(1:ndir), xp(3), told, tnew
217 
218  if (const_dt_particles > 0) then
219  dt_p = const_dt_particles
220  return
221  end if
222 
223  dt_p = dtsave_particles
224 
225  ! Make sure the user-defined particle movement doesn't break communication
226  if (associated(usr_particle_position)) then
227  xp = partp%self%x
228  told = partp%self%time
229  tnew = told+dt_p
230  call usr_particle_position(xp, partp%self%index, told, tnew)
231  do while (.not. point_in_igrid_ghostc(xp,partp%igrid,1))
232  dt_p = dt_p/10.d0
233  xp = partp%self%x
234  tnew = told+dt_p
235  call usr_particle_position(xp, partp%self%index, told, tnew)
236  end do
237  end if
238 
239  ! Make sure we don't advance beyond end_time
240  call limit_dt_endtime(end_time - partp%self%time, dt_p)
241 
242  end function sample_get_particle_dt
243 
244 end module mod_particle_sample
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
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.
logical saveprim
If true, convert from conservative to primitive variables in output.
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.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision, dimension(ndim) dxlevel
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.
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 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.
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
procedure(sub_integrate), pointer particles_integrate
procedure(sub_fill_gridvars), pointer particles_fill_gridvars
integer ndefpayload
Number of default payload variables for a particle.
logical function point_in_igrid_ghostc(x, igrid, ngh)
Quick check if particle coordinate is inside igrid (ghost cells included, except the last ngh)
integer igrid_working
set the current igrid for the particle integrator:
subroutine interpolate_var(igrid, ixIL, ixOL, gf, x, xloc, gfloc)
Scattered sampling based on fixed- or moving-particle interpolation By Fabio Bacchini (2020)
double precision function sample_get_particle_dt(partp, end_time)
subroutine, public sample_init()
subroutine sample_integrate_particles(end_time)
subroutine sample_update_payload(igrid, xpart, upart, qpart, mpart, mypayload, mynpayload, particle_time)
Payload update.
subroutine, public sample_create_particles()
Module with all the methods that users can customize in AMRVAC.
procedure(particle_position), pointer usr_particle_position
procedure(check_particle), pointer usr_check_particle
procedure(create_particles), pointer usr_create_particles
procedure(update_payload), pointer usr_update_payload