Afivo 0.3
All Classes Namespaces Functions Variables Pages
advection.f90
1#include "../src/cpp_macros.h"
2!> \example advection_v2.f90
3!>
4!> An advection example
5program advection
6 use m_af_all
7
8 implicit none
9
10 integer, parameter :: box_size = 8
11 integer :: i_phi
12 integer :: i_err
13 integer :: i_flux
14 integer, parameter :: sol_type = 2
15 real(dp), parameter :: domain_len = 2 * acos(-1.0_dp)
16
17 type(af_t) :: tree
18 type(ref_info_t) :: refine_info
19 integer :: refine_steps, output_cnt
20 real(dp) :: dt, dt_lim, time, time_prev_refine
21 real(dp) :: end_time, err, sum_err2
22 real(dp) :: sum_phi, sum_phi_t0
23 real(dp) :: dt_amr, dt_output
24 real(dp) :: velocity(NDIM)
25 real(dp) :: cfl_number
26 integer :: integrator
27 character(len=100) :: fname
28 character(len=20) :: flux_method
29
30 print *, "Running advection_" // dimname // ""
31 print *, "Number of threads", af_get_max_threads()
32
33 flux_method = "generic"
34
35 if (command_argument_count() > 0) then
36 call get_command_argument(1, flux_method)
37 end if
38
39 integrator = af_heuns_method
40
41 call af_add_cc_variable(tree, "phi", ix=i_phi, &
42 n_copies=af_advance_num_steps(integrator))
43 call af_add_cc_variable(tree, "err", ix=i_err)
44 call af_add_fc_variable(tree, "flux", ix=i_flux)
45
46 call af_set_cc_methods(tree, i_phi, af_bc_neumann_zero, &
47 prolong=af_prolong_limit)
48
49 ! Initialize tree
50 call af_init(tree, & ! Tree to initialize
51 box_size, & ! A box contains box_size**DIM cells
52 [dtimes(domain_len)], &
53 [dtimes(box_size)], &
54 periodic=[dtimes(.true.)])
55
56 output_cnt = 0
57 time = 0
58 dt_amr = 0.01_dp
59 dt_output = 0.5_dp
60 end_time = 5.0_dp
61 velocity(:) = -0.5_dp
62 velocity(1) = 1.0_dp
63 cfl_number = 0.5_dp
64
65 ! Set up the initial conditions
66 refine_steps=0
67
68 do
69 refine_steps=refine_steps+1
70 ! Set initial conditions on all boxes
71 call af_loop_box(tree, set_initial_condition)
72
73 ! Fill ghost cells for variables i_phi
74 call af_gc_tree(tree, [i_phi])
75
76 ! Adjust the refinement of a tree using refine_routine
77 call af_adjust_refinement(tree, refine_routine, refine_info, 1)
78
79 ! If no new boxes have been added, exit the loop
80 if (refine_info%n_add == 0) exit
81 end do
82
83 call af_print_info(tree)
84
85 ! Restrict the initial conditions
86 call af_restrict_tree(tree, [i_phi])
87
88 ! Fill ghost cells for variables i_phi on the sides of all boxes
89 call af_gc_tree(tree, [i_phi])
90
91 call af_tree_sum_cc(tree, i_phi, sum_phi_t0)
92
93 dt = cfl_number / (sum(abs(velocity/af_min_dr(tree))) + epsilon(1.0_dp))
94 time_prev_refine = time
95
96 ! Starting simulation
97 do
98 if (output_cnt * dt_output <= time) then
99 output_cnt = output_cnt + 1
100 write(fname, "(A,I0)") "output/advection_v2_" // dimname // "_", output_cnt
101
102 ! Call procedure set_error (see below) for each box in tree, with argument time
103 call af_loop_box_arg(tree, set_error, [time])
104
105 ! Write the cell centered data of tree to a vtk unstructured file fname.
106 ! Only the leaves of the tree are used
107 call af_write_silo(tree, trim(fname), output_cnt, time)
108
109 ! Find maximum and minimum values of cc(..., i_err) and cc(..., i_phi).
110 ! By default, only loop over leaves, and ghost cells are not used.
111 call af_tree_maxabs_cc(tree, i_err, err)
112 call af_tree_sum_cc(tree, i_phi, sum_phi)
113 call af_tree_sum_cc(tree, i_err, sum_err2, power=2)
114 write(*,"(3(A,Es16.8))") &
115 " max error:", err, &
116 " mean error:", sqrt(sum_err2/af_total_volume(tree)), &
117 " conservation error: ", (sum_phi - sum_phi_t0) / sum_phi_t0
118 end if
119
120 if (time > end_time) exit
121
122 call af_advance(tree, dt, dt_lim, time, [i_phi], integrator, forward_euler)
123 dt = cfl_number * dt_lim
124
125 if (time > time_prev_refine + dt_amr) then
126 call af_restrict_tree(tree, [i_phi])
127 call af_gc_tree(tree, [i_phi])
128 call af_adjust_refinement(tree, refine_routine, refine_info, 1)
129 time_prev_refine = time
130 end if
131 end do
132
133contains
134
135 !> [refine_routine]
136 !> Set refinement flags for box
137 subroutine refine_routine(box, cell_flags)
138 type(box_t), intent(in) :: box
139 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
140 real(dp) :: diff
141 integer :: IJK, nc
142
143 nc = box%n_cell
144
145 do kji_do(1,nc)
146#if NDIM == 1
147 diff = abs(box%dr(1) * (box%cc(i+1, i_phi) + &
148 box%cc(i-1, i_phi) - 2 * box%cc(i, i_phi)))
149#elif NDIM == 2
150 diff = abs(box%dr(1) * (box%cc(i+1, j, i_phi) + &
151 box%cc(i-1, j, i_phi) - 2 * box%cc(i, j, i_phi)) + &
152 box%dr(2) * (box%cc(i, j+1, i_phi) + &
153 box%cc(i, j-1, i_phi) - 2 * box%cc(i, j, i_phi)))
154#elif NDIM == 3
155 diff = abs(box%dr(1) * (box%cc(i+1, j, k, i_phi) + &
156 box%cc(i-1, j, k, i_phi) - 2 * box%cc(i, j, k, i_phi)) + &
157 box%dr(2) * (box%cc(i, j+1, k, i_phi) + &
158 box%cc(i, j-1, k, i_phi) - 2 * box%cc(i, j, k, i_phi)) + &
159 box%dr(3) * (box%cc(i, j, k+1, i_phi) + &
160 box%cc(i, j, k-1, i_phi) - 2 * box%cc(i, j, k, i_phi)))
161#endif
162
163 if (box%lvl < 2 .or. diff > 2.0e-3_dp .and. box%lvl < 5) then
164 cell_flags(ijk) = af_do_ref
165 else if (diff < 0.1_dp * 2.0e-3_dp) then
166 cell_flags(ijk) = af_rm_ref
167 else
168 cell_flags(ijk) = af_keep_ref
169 end if
170 end do; close_do
171 end subroutine refine_routine
172 !> [refine_routine]
173
174 !> This routine sets the initial conditions for each box
175 subroutine set_initial_condition(box)
176 type(box_t), intent(inout) :: box
177 integer :: IJK, nc
178 real(dp) :: rr(NDIM)
179
180 nc = box%n_cell
181 do kji_do(0,nc+1)
182 rr = af_r_cc(box, [ijk])
183 box%cc(ijk, i_phi) = solution(rr, 0.0_dp)
184 end do; close_do
185 end subroutine set_initial_condition
186
187 !> This routine computes the error in i_phi
188 subroutine set_error(box, time)
189 type(box_t), intent(inout) :: box
190 real(dp), intent(in) :: time(:)
191 integer :: IJK, nc
192 real(dp) :: rr(NDIM)
193
194 nc = box%n_cell
195 do kji_do(1,nc)
196 rr = af_r_cc(box, [ijk])
197 box%cc(ijk, i_err) = &
198 box%cc(ijk, i_phi) - solution(rr, time(1))
199 end do; close_do
200 end subroutine set_error
201
202 !> This routine calculates the analytic solution in point rr
203 function solution(rr, t) result(sol)
204 real(dp), intent(in) :: rr(NDIM), t
205 real(dp) :: sol, rr_t(NDIM)
206
207 rr_t = rr - velocity * t
208
209 select case (sol_type)
210 case (1)
211#if NDIM > 1
212 sol = sin(0.5_dp * rr_t(1))**8 * cos(0.5_dp * rr_t(2))**8
213#else
214 sol = sin(0.5_dp * rr_t(1))**8
215#endif
216 case (2)
217 rr_t = modulo(rr_t, domain_len) / domain_len
218 if (norm2(rr_t - 0.5_dp) < 0.1_dp) then
219 sol = 1.0_dp
220 else
221 sol = 0.0_dp
222 end if
223 end select
224 end function solution
225
226 subroutine forward_euler(tree, dt, dt_stiff, dt_lim, time, s_deriv, n_prev, &
227 s_prev, w_prev, s_out, i_step, n_steps)
228 type(af_t), intent(inout) :: tree
229 real(dp), intent(in) :: dt
230 real(dp), intent(in) :: dt_stiff !< Time step for stiff terms
231 real(dp), intent(in) :: time
232 real(dp), intent(inout) :: dt_lim
233 integer, intent(in) :: s_deriv
234 integer, intent(in) :: n_prev !< Number of previous states
235 integer, intent(in) :: s_prev(n_prev) !< Previous states
236 real(dp), intent(in) :: w_prev(n_prev) !< Weights of previous states
237 integer, intent(in) :: s_out
238 integer, intent(in) :: i_step, n_steps
239 integer :: lvl, i, id
240 real(dp) :: all_dt(1)
241
242 select case (flux_method)
243 case ("generic")
244 call flux_generic_tree(tree, 1, [i_phi], s_deriv, [i_flux], dt_lim, &
245 max_wavespeed, get_flux, flux_dummy_modify, flux_dummy_line_modify, &
246 flux_dummy_conversion, flux_dummy_conversion, af_limiter_koren_t)
247 case ("upwind")
248 call flux_upwind_tree(tree, 1, [i_phi], s_deriv, [i_flux], &
249 1, all_dt, flux_upwind, flux_direction, &
250 flux_dummy_line_modify, af_limiter_koren_t)
251 dt_lim = all_dt(1)
252 case ("custom")
253 ! Ensure ghost cells near refinement boundaries can be properly filled
254 call af_restrict_ref_boundary(tree, [i_phi+s_deriv])
255
256 ! Compute fluxes
257 !$omp parallel private(lvl, i, id)
258 do lvl = 1, tree%highest_lvl
259 !$omp do
260 do i = 1, size(tree%lvls(lvl)%leaves)
261 id = tree%lvls(lvl)%leaves(i)
262 call fluxes_koren(tree, id, s_deriv)
263 end do
264 !$omp end do
265 end do
266 !$omp end parallel
267
268 ! Restrict fluxes from children to parents on refinement boundaries.
269 call af_consistent_fluxes(tree, [i_flux])
270
271 dt_lim = 1/(sum(abs(velocity/af_min_dr(tree))) + epsilon(1.0_dp))
272 case default
273 error stop "Unknown flux_method, choices: generic, upwind, custom"
274 end select
275
276 call flux_update_densities(tree, dt, 1, [i_phi], 1, [i_phi], [i_flux], &
277 s_deriv, n_prev, s_prev, w_prev, s_out, flux_dummy_source, 0, all_dt)
278
279 end subroutine forward_euler
280
281 !> This routine computes the x-fluxes and y-fluxes interior (advective part)
282 !> with the Koren limiter
283 subroutine fluxes_koren(tree, id, s_deriv)
285 type(af_t), intent(inout) :: tree
286 integer, intent(in) :: id
287 integer, intent(in) :: s_deriv
288 integer :: nc, idim
289 real(dp), allocatable :: cc(DTIMES(:), :)
290 real(dp), allocatable :: v(DTIMES(:), :)
291
292 nc = tree%boxes(id)%n_cell
293 allocate(cc(dtimes(-1:nc+2), 1))
294 allocate(v(dtimes(1:nc+1), ndim))
295
296 call af_gc2_box(tree, id, [i_phi+s_deriv], cc)
297
298 do idim = 1, ndim
299 v(dtimes(:), idim) = velocity(idim)
300 end do
301
302#if NDIM == 1
303 call flux_koren_1d(cc(dtimes(:), 1), v, nc, 2)
304#elif NDIM == 2
305 call flux_koren_2d(cc(dtimes(:), 1), v, nc, 2)
306#elif NDIM == 3
307 call flux_koren_3d(cc(dtimes(:), 1), v, nc, 2)
308#endif
309
310 tree%boxes(id)%fc(dtimes(:), :, i_flux) = v
311
312 end subroutine fluxes_koren
313
314 subroutine max_wavespeed(n_values, n_var, flux_dim, u, w)
315 integer, intent(in) :: n_values !< Number of cell faces
316 integer, intent(in) :: n_var !< Number of variables
317 integer, intent(in) :: flux_dim !< In which dimension fluxes are computed
318 real(dp), intent(in) :: u(n_values, n_var) !< Primitive variables
319 real(dp), intent(out) :: w(n_values) !< Maximum speed
320
321 w = abs(velocity(flux_dim))
322 end subroutine max_wavespeed
323
324 subroutine get_flux(n_values, n_var, flux_dim, u, flux, box, line_ix, s_deriv)
325 integer, intent(in) :: n_values !< Number of cell faces
326 integer, intent(in) :: n_var !< Number of variables
327 integer, intent(in) :: flux_dim !< In which dimension fluxes are computed
328 real(dp), intent(in) :: u(n_values, n_var)
329 real(dp), intent(out) :: flux(n_values, n_var)
330 type(box_t), intent(in) :: box
331 integer, intent(in) :: line_ix(NDIM-1)
332 integer, intent(in) :: s_deriv
333
334 flux = u * velocity(flux_dim)
335 end subroutine get_flux
336
337 subroutine flux_upwind(nf, n_var, flux_dim, u, flux, cfl_sum, &
338 n_other_dt, other_dt, box, line_ix, s_deriv)
339 integer, intent(in) :: nf !< Number of cell faces
340 integer, intent(in) :: n_var !< Number of variables
341 integer, intent(in) :: flux_dim !< In which dimension fluxes are computed
342 real(dp), intent(in) :: u(nf, n_var) !< Face values
343 real(dp), intent(out) :: flux(nf, n_var) !< Computed fluxes
344 !> Terms per cell-center to be added to CFL sum, see flux_upwind_box
345 real(dp), intent(out) :: cfl_sum(nf-1)
346 integer, intent(in) :: n_other_dt !< Number of non-cfl time step restrictions
347 real(dp), intent(inout) :: other_dt(n_other_dt) !< Non-cfl time step restrictions
348 type(box_t), intent(in) :: box !< Current box
349 integer, intent(in) :: line_ix(NDIM-1) !< Index of line for dim /= flux_dim
350 integer, intent(in) :: s_deriv !< State to compute derivatives from
351 real(dp) :: tmp
352
353 flux = u * velocity(flux_dim)
354
355 tmp = abs(velocity(flux_dim)) / box%dr(flux_dim)
356 cfl_sum = tmp
357 end subroutine flux_upwind
358
359 subroutine flux_direction(box, line_ix, s_deriv, n_var, flux_dim, direction_positive)
360 type(box_t), intent(in) :: box !< Current box
361 integer, intent(in) :: line_ix(NDIM-1) !< Index of line for dim /= flux_dim
362 integer, intent(in) :: s_deriv !< State to compute derivatives from
363 integer, intent(in) :: n_var !< Number of variables
364 integer, intent(in) :: flux_dim !< In which dimension fluxes are computed
365 !> True means positive flow (to the "right"), false to the left
366 logical, intent(out) :: direction_positive(box%n_cell+1, n_var)
367
368 direction_positive(:, 1) = (velocity(flux_dim) > 0)
369 end subroutine flux_direction
370
371end program
Module which contains all Afivo modules, so that a user does not have to include them separately.
Definition: m_af_all.f90:3
Module containing a couple flux schemes for solving hyperbolic problems explicitly,...