Toy example showing how to simulate gravitating particles.
1#include "../src/cpp_macros.h"
2
3
4
5program particles_gravity
7
8 implicit none
9
10 integer :: n
11 integer, parameter :: box_size = 8
12 integer, parameter :: n_particles = 100*1000
13 integer, parameter :: particles_per_cell = 100
14 integer, parameter :: max_refinement_lvl = 5
15 integer :: i_phi
16 integer :: i_rho
17 integer :: i_tmp
18 integer :: i_f(NDIM)
19
20 real(dp), parameter :: force_to_accel = 1.0_dp
21
22 real(dp), parameter :: domain_len = 1.0_dp
23 real(dp), parameter :: mean_density = 1.0_dp
24 real(dp), parameter :: particle_weight = domain_len**ndim / n_particles
25 real(dp), parameter :: t_end = 0.1_dp
26 real(dp), parameter :: dt_output = 5e-2_dp
27 real(dp) :: dt = 1.0e-2_dp
28 real(dp) :: dt_prev
29
30 real(dp), allocatable :: coordinates(:, :)
31 real(dp), allocatable :: velocities(:, :)
32 integer, allocatable :: ids(:)
33 real(dp), allocatable :: weights(:)
34
35 integer :: n_output
36 type(af_t) :: tree
37 type(mg_t) :: mg
38 type(ref_info_t) :: refine_info
39 integer :: count_rate, wc_start, wc_end
40 character(len=100) :: fname
41 real(dp) :: max_vel, time
42
43 print *, "Running particles_gravity_" // dimname // ""
44 print *, "Number of threads", af_get_max_threads()
45
46 call af_add_cc_variable(tree, "phi", ix=i_phi)
47 call af_add_cc_variable(tree, "rho", ix=i_rho)
48 call af_add_cc_variable(tree, "tmp", ix=i_tmp)
49 call af_add_cc_variable(tree, "fx", ix=i_f(1))
50#if NDIM > 1
51 call af_add_cc_variable(tree, "fy", ix=i_f(2))
52#endif
53#if NDIM == 3
54 call af_add_cc_variable(tree, "fz", ix=i_f(3))
55#endif
56
57
58 call af_init(tree, &
59 box_size, &
60 [dtimes(domain_len)], &
61 [dtimes(box_size)], &
62 periodic=[dtimes(.true.)])
63
64 mg%i_phi = i_phi
65 mg%i_rhs = i_rho
66 mg%i_tmp = i_tmp
67 mg%sides_bc => af_bc_neumann_zero
68 mg%subtract_mean = .true.
69
70 call mg_init(tree, mg)
71
72
73
74 call af_set_cc_methods(tree, i_rho, af_bc_neumann_zero, af_gc_interp)
75 do n = 1, ndim
76 call af_set_cc_methods(tree, i_f(n), af_bc_neumann_zero, af_gc_interp)
77 end do
78
79 allocate(coordinates(ndim, n_particles))
80 allocate(velocities(ndim, n_particles))
81 allocate(weights(n_particles))
82 allocate(ids(n_particles))
83
84 call random_number(coordinates(:, :))
85 call random_number(velocities(:, :))
86 velocities(:, :) = 0.0_dp
87 weights(:) = particle_weight
88
89 do n = 1, 100
90 call af_tree_clear_cc(tree, i_rho)
91 call af_particles_to_grid(tree, i_rho, n_particles, get_id, get_rw, 1)
92 call af_adjust_refinement(tree, refine_routine, refine_info, 2)
93 if (refine_info%n_add == 0) exit
94 end do
95
96
97 call mg_fas_fmg(tree, mg, .false., .false.)
98 call mg_fas_fmg(tree, mg, .false., .true.)
99
100 time = 0.0_dp
101 n = 0
102 n_output = 0
103 dt_prev = dt
104
105 call system_clock(wc_start, count_rate)
106 do while (time < t_end)
107 n = n + 1
108 call push_particles(coordinates, velocities, n_particles, dt)
109
110
111 call af_tree_clear_cc(tree, i_rho)
112 call af_particles_to_grid(tree, i_rho, n_particles, get_id, get_rw, 1)
113 call af_loop_box(tree, subtract_mean_density)
114
115 call mg_fas_vcycle(tree, mg, .false.)
116 call mg_fas_vcycle(tree, mg, .true.)
117
118 call af_loop_box(tree, compute_forces)
119 call af_gc_tree(tree, i_f)
120
121 call update_velocities(tree, coordinates, velocities, &
122 ids, n_particles, dt)
123
124 time = time + dt
125 max_vel = sqrt(maxval(sum(velocities**2, dim=1)))
126 dt = min(1.1 * dt_prev, dt_output, &
127 0.5_dp * af_min_dr(tree) / max_vel)
128 dt_prev = dt
129 call compute_total_energy(tree, coordinates, velocities, &
130 ids, n_particles)
131
132 if (modulo(n, 5) == 0) then
133 call af_adjust_refinement(tree, refine_routine, refine_info, 2)
134 end if
135
136 if (n_output * dt_output <= time) then
137 call af_tree_clear_cc(tree, i_rho)
138 call af_particles_to_grid(tree, i_rho, n_particles, get_id, get_rw, 1)
139 n_output = n_output + 1
140 write(fname, "(A,I4.4)") "output/particles_gravity_" // dimname // "_", n_output
141 call af_write_silo(tree, fname, n_output, time)
142 print *, "Time", time, "n_cells", af_num_cells_used(tree)
143 end if
144 end do
145
146 call system_clock(wc_end, count_rate)
147 print *, "Total time: ", (wc_end-wc_start) / real(count_rate, dp), " seconds"
148
149contains
150
151 subroutine compute_total_energy(tree, x, v, ids, n_particles)
152 type(af_t), intent(in) :: tree
153 integer, intent(in) :: n_particles
154 real(dp), intent(in) :: x(NDIM, n_particles)
155 real(dp), intent(in) :: v(NDIM, n_particles)
156 integer, intent(inout) :: ids(n_particles)
157 real(dp) :: sum_ken, sum_phi, phi(1)
158 integer :: n
159 logical :: success
160
161 sum_ken = 0.0_dp
162 sum_phi = 0.0_dp
163
164 do n = 1, n_particles
165 phi = af_interp1(tree, x(:, n), [i_phi], success, ids(n))
166 if (.not. success) error stop "compute_total_energy: interpolation error"
167 sum_ken = sum_ken + 0.5_dp * sum(v(:, n)**2)
168 sum_phi = sum_phi + 0.5_dp * phi(1)
169 end do
170 print *, "sum energy", sum_ken, sum_phi, sum_ken+sum_phi
171 end subroutine compute_total_energy
172
173 subroutine push_particles(x, v, n_particles, dt)
174 integer, intent(in) :: n_particles
175 real(dp), intent(inout) :: x(NDIM, n_particles)
176 real(dp), intent(in) :: v(NDIM, n_particles)
177 real(dp), intent(in) :: dt
178 integer :: n
179
180
181 do n = 1, n_particles
182 x(:, n) = x(:, n) + dt * v(:, n)
183 x(:, n) = modulo(x(:, n), domain_len)
184 end do
185
186 end subroutine push_particles
187
188 subroutine update_velocities(tree, x, v, ids, n_particles, dt)
189 type(af_t), intent(in) :: tree
190 integer, intent(in) :: n_particles
191 real(dp), intent(in) :: x(NDIM, n_particles)
192 real(dp), intent(inout) :: v(NDIM, n_particles)
193 integer, intent(inout) :: ids(:)
194 real(dp), intent(in) :: dt
195 integer :: n
196 real(dp) :: a(NDIM)
197 logical :: success
198
199
200 do n = 1, n_particles
201 a = af_interp1(tree, x(:, n), i_f, success, id_guess=ids(n))
202 if (.not. success) error stop "update velocities: interpolation error"
203 v(:, n) = v(:, n) + dt * a * force_to_accel
204 end do
205
206 end subroutine update_velocities
207
208 subroutine refine_routine(box, cell_flags)
209 type(box_t), intent(in) :: box
210 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
211 real(dp) :: rho_refine, rho_derefine
212 integer :: nc
213
214 nc = box%n_cell
215 rho_refine = particles_per_cell * particle_weight / product(box%dr)
216 rho_derefine = rho_refine / (2**ndim * 4)
217
218 where (box%cc(dtimes(1:nc), i_rho) > rho_refine)
219 cell_flags(dtimes(:)) = af_do_ref
220 elsewhere (box%cc(dtimes(1:nc), i_rho) < rho_derefine)
221 cell_flags(dtimes(:)) = af_rm_ref
222 elsewhere
223 cell_flags(dtimes(:)) = af_keep_ref
224 end where
225
226 if (box%lvl >= max_refinement_lvl) then
227 where (cell_flags == af_do_ref)
228 cell_flags = af_keep_ref
229 end where
230 end if
231
232 end subroutine refine_routine
233
234 subroutine compute_forces(box)
235 type(box_t), intent(inout) :: box
236 integer :: IJK, nc
237 real(dp) :: fac(NDIM)
238
239 nc = box%n_cell
240 fac = -0.5_dp / box%dr
241
242 do kji_do(1, nc)
243#if NDIM == 1
244 box%cc(i, i_f(1)) = fac(1) * &
245 (box%cc(i+1, i_phi) - box%cc(i-1, i_phi))
246#elif NDIM == 2
247 box%cc(i, j, i_f(1)) = fac(1) * &
248 (box%cc(i+1, j, i_phi) - box%cc(i-1, j, i_phi))
249 box%cc(i, j, i_f(2)) = fac(2) * &
250 (box%cc(i, j+1, i_phi) - box%cc(i, j-1, i_phi))
251#elif NDIM == 3
252 box%cc(i, j, k, i_f(1)) = fac(1) * &
253 (box%cc(i+1, j, k, i_phi) - box%cc(i-1, j, k, i_phi))
254 box%cc(i, j, k, i_f(2)) = fac(2) * &
255 (box%cc(i, j+1, k, i_phi) - box%cc(i, j-1, k, i_phi))
256 box%cc(i, j, k, i_f(3)) = fac(3) * &
257 (box%cc(i, j, k+1, i_phi) - box%cc(i, j, k-1, i_phi))
258#endif
259 end do; close_do
260
261 end subroutine compute_forces
262
263 subroutine subtract_mean_density(box)
264 type(box_t), intent(inout) :: box
265
266 box%cc(dtimes(:), i_rho) = box%cc(dtimes(:), i_rho) - mean_density
267 end subroutine subtract_mean_density
268
269 subroutine get_id(n, id)
270 integer, intent(in) :: n
271 integer, intent(out) :: id
272
273 id = af_get_id_at(tree, coordinates(:, n), guess=ids(n))
274 ids(n) = id
275 end subroutine get_id
276
277 subroutine get_rw(n, r, w)
278 integer, intent(in) :: n
279 real(dp), intent(out) :: r(NDIM)
280 real(dp), intent(out) :: w
281
282 r = coordinates(:, n)
283 w = weights(n)
284 end subroutine get_rw
285
286end program particles_gravity
Module which contains all Afivo modules, so that a user does not have to include them separately.