Example showing how to map particles to a grid.
1#include "../src/cpp_macros.h"
2
3
4
5program particles_to_grid
7
8 implicit none
9
10 integer, parameter :: box_size = 8
11 integer :: i_phi, i_tmp
12 integer, parameter :: n_particles = 1000*1000
13 real(dp), parameter :: r_max(NDIM) = 1.0_dp
14 real(dp), allocatable :: coordinates(:, :), weights(:)
15 integer :: coordinate_system = af_xyz
16 logical :: cylindrical = .false.
17 logical :: use_tmp_var = .false.
18
19 type(af_t) :: tree
20 type(ref_info_t) :: refine_info
21 integer :: n
22 integer :: count_rate, t_start, t_end
23 real(dp) :: sum_density
24 character(len=10) :: tmp_string
25
26 if (command_argument_count() > 0) then
27 call get_command_argument(1, tmp_string)
28 read(tmp_string, *) cylindrical
29 end if
30
31 if (command_argument_count() > 1) then
32 call get_command_argument(2, tmp_string)
33 read(tmp_string, *) use_tmp_var
34 end if
35
36 if (cylindrical) coordinate_system = af_cyl
37
38 print *, "Running particles_to_grid_" // dimname // ""
39 print *, "Number of threads", af_get_max_threads()
40 print *, "Cylindrical", cylindrical
41 print *, "Number of particles", n_particles
42
43 call af_add_cc_variable(tree, "phi", ix=i_phi)
44 call af_add_cc_variable(tree, "tmp", ix=i_tmp)
45
46 call af_init(tree, box_size, r_max, [dtimes(box_size)], &
47 coord=coordinate_system)
48
49 call af_set_cc_methods(tree, i_phi, af_bc_neumann_zero, af_gc_interp)
50
51 do
52 call af_adjust_refinement(tree, refine_routine, refine_info, 0)
53 if (refine_info%n_add == 0) exit
54 end do
55
56 allocate(coordinates(ndim, n_particles))
57 allocate(weights(n_particles))
58
59 call random_number(coordinates(:, :))
60 weights(:) = 1.0_dp / n_particles
61
62 if (cylindrical) then
63 coordinates(1, :) = sqrt(coordinates(1, :))
64 end if
65
66 do n = 1, n_particles
67 coordinates(:, n) = 0.6_dp * coordinates(:, n) * r_max
68 end do
69
70 call system_clock(t_start, count_rate)
71
72 if (use_tmp_var) then
73 call af_particles_to_grid(tree, i_phi, n_particles, &
74 get_id, get_rw, 0, iv_tmp=i_tmp)
75 else
76 call af_particles_to_grid(tree, i_phi, n_particles, &
77 get_id, get_rw, 0)
78 end if
79
80 call system_clock(t_end, count_rate)
81 call af_write_silo(tree, "output/particles_to_grid_" // dimname // "_0", 1, 0.0_dp)
82 call af_tree_sum_cc(tree, i_phi, sum_density)
83 print *, "zeroth order: ", (t_end-t_start) / real(count_rate, dp), " seconds"
84 print *, "integrated density: ", sum_density
85
86 call af_tree_clear_cc(tree, i_phi)
87
88 call system_clock(t_start, count_rate)
89
90 if (use_tmp_var) then
91 call af_particles_to_grid(tree, i_phi, n_particles, &
92 get_id, get_rw, 1, iv_tmp=i_tmp)
93 else
94 call af_particles_to_grid(tree, i_phi, n_particles, &
95 get_id, get_rw, 1)
96 end if
97
98 call system_clock(t_end, count_rate)
99 call af_write_silo(tree, "output/particles_to_grid_" // dimname // "_1", 1, 0.0_dp)
100 call af_tree_sum_cc(tree, i_phi, sum_density)
101 print *, "first order: ", (t_end-t_start) / real(count_rate, dp), " seconds"
102 print *, "integrated density: ", sum_density
103
104contains
105
106 subroutine refine_routine(box, cell_flags)
107 type(box_t), intent(in) :: box
108 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
109
110 if (all(box%r_min < 0.5_dp * r_max) .and. box%lvl <= 5) then
111 cell_flags(dtimes(:)) = af_do_ref
112 else
113 cell_flags(dtimes(:)) = af_keep_ref
114 end if
115 end subroutine refine_routine
116
117 subroutine get_id(n, id)
118 integer, intent(in) :: n
119 integer, intent(out) :: id
120
121 id = af_get_id_at(tree, coordinates(:, n))
122 end subroutine get_id
123
124 subroutine get_rw(n, r, w)
125 integer, intent(in) :: n
126 real(dp), intent(out) :: r(NDIM)
127 real(dp), intent(out) :: w
128
129 r = coordinates(:, n)
130 w = weights(n)
131 end subroutine get_rw
132
133end program particles_to_grid
Module which contains all Afivo modules, so that a user does not have to include them separately.
