Afivo 0.3
All Classes Namespaces Functions Variables Pages
m_af_core.f90
1!> This module contains the core routines of Afivo, namely those that deal with
2!> initializing and changing the quadtree/octree mesh.
4#include "cpp_macros.h"
5 use m_af_types
6
7 implicit none
8 private
9
10 public :: af_add_cc_variable
11 public :: af_add_fc_variable
12 public :: af_find_cc_variable
13 public :: af_find_fc_variable
14 public :: af_init
15 public :: af_set_cc_methods
16 public :: af_init_box
17 public :: af_destroy
18 public :: af_adjust_refinement
19 public :: af_refine_up_to_lvl
20 public :: af_consistent_fluxes
21
22contains
23
24 !> Add cell-centered variable
25 !> @todo ix as third argument?
26 subroutine af_add_cc_variable(tree, name, write_out, n_copies, &
27 ix, write_binary)
28 !> Tree to add variable to
29 type(af_t), intent(inout) :: tree
30 !> Name of the variable
31 character(len=*), intent(in) :: name
32 !> Include variable in output
33 logical, intent(in), optional :: write_out
34 !> Include variable in binary output (for restarting)
35 logical, intent(in), optional :: write_binary
36 !> How many copies of variable to store (default: 1)
37 integer, intent(in), optional :: n_copies
38 !> On output: index of variable
39 integer, intent(out), optional :: ix
40
41 integer :: n, ncpy
42 logical :: writeout, writebin
43
44 ncpy = 1; if (present(n_copies)) ncpy = n_copies
45 writeout = .true.; if (present(write_out)) writeout = write_out
46 writebin = .true.; if (present(write_binary)) writebin = write_binary
47
48 if (ncpy < 1) error stop "af_add_cc_variable: n_copies < 1"
49
50 if (tree%n_var_cell + ncpy > af_max_num_vars) then
51 print *, "af_max_num_vars:", af_max_num_vars
52 print *, "Cannot add ", name
53 error stop "Too many cc variables"
54 end if
55
56 do n = 1, ncpy
57 tree%n_var_cell = tree%n_var_cell + 1
58 if (n == 1) then
59 if (present(ix)) ix = tree%n_var_cell
60 tree%cc_names(tree%n_var_cell) = name
61 tree%cc_write_output(tree%n_var_cell) = writeout
62 tree%cc_write_binary(tree%n_var_cell) = writebin
63 tree%cc_num_copies(tree%n_var_cell) = ncpy
64 else
65 write(tree%cc_names(tree%n_var_cell), "(A,I0)") &
66 trim(name) // '_', n
67 tree%cc_write_output(tree%n_var_cell) = .false.
68 tree%cc_write_binary(tree%n_var_cell) = .false.
69 tree%cc_num_copies(tree%n_var_cell) = 0
70 end if
71 end do
72
73 end subroutine af_add_cc_variable
74
75 !> Add face-centered variable
76 subroutine af_add_fc_variable(tree, name, ix, write_binary)
77 !> Tree to add variable to
78 type(af_t), intent(inout) :: tree
79 !> Name of the variable
80 character(len=*), intent(in) :: name
81 !> On output: index of variable
82 integer, intent(out), optional :: ix
83 !> Include variable in binary output
84 logical, intent(in), optional :: write_binary
85 logical :: writebin
86
87 writebin = .true.; if (present(write_binary)) writebin = write_binary
88
89 if (tree%n_var_face + 1 > af_max_num_vars) then
90 print *, "af_max_num_vars:", af_max_num_vars
91 print *, "Cannot add ", name
92 error stop "Too many fc variables"
93 end if
94
95 tree%n_var_face = tree%n_var_face + 1
96 tree%fc_names(tree%n_var_face) = name
97 tree%fc_write_binary(tree%n_var_face) = writebin
98 if (present(ix)) ix = tree%n_var_face
99 end subroutine af_add_fc_variable
100
101 !> Find index of cell-centered variable
102 integer function af_find_cc_variable(tree, name)
103 type(af_t), intent(in) :: tree
104 character(len=*), intent(in) :: name
105 integer :: n
106
107 do n = 1, tree%n_var_cell
108 if (tree%cc_names(n) == name) exit
109 end do
110
111 if (n == tree%n_var_cell+1) then
112 print *, "variable name: ", trim(name)
113 error stop "af_find_cc_variable: variable not found"
114 end if
115
116 af_find_cc_variable = n
117 end function af_find_cc_variable
118
119 !> Find index of face-centered variable
120 integer function af_find_fc_variable(tree, name)
121 type(af_t), intent(in) :: tree
122 character(len=*), intent(in) :: name
123 integer :: n
124
125 do n = 1, tree%n_var_face
126 if (tree%fc_names(n) == name) exit
127 end do
128
129 if (n == tree%n_var_face+1) then
130 print *, "variable name: ", trim(name)
131 error stop "af_find_fc_variable: variable not found"
132 end if
133
134 af_find_fc_variable = n
135 end function af_find_fc_variable
136
137 !> Initialize a NDIM-d octree/quadtree grid
138 subroutine af_init(tree, n_cell, r_max, grid_size, periodic, r_min, coord, &
139 mem_limit_gb, box_limit)
140 type(af_t), intent(inout) :: tree !< The tree to initialize
141 integer, intent(in) :: n_cell !< Boxes have n_cell^dim cells
142 real(dp), intent(in) :: r_max(ndim) !< Maximal coordinates of the domain
143 integer, intent(in) :: grid_size(ndim) !< Size of the coarse grid
144 logical, intent(in), optional :: periodic(ndim) !< True for periodic dimensions
145 real(dp), intent(in), optional :: r_min(ndim) !< Lowest coordinate, default is (0., 0., 0.)
146 integer, intent(in), optional :: coord !< Select coordinate type
147 real(dp), intent(in), optional :: mem_limit_gb !< Memory limit in GByte
148 !> Maximum number of boxes (overrides mem_limit_gb)
149 integer, intent(in), optional :: box_limit
150
151 real(dp) :: r_min_a(ndim), gb_limit
152 integer :: lvl, coord_a, box_bytes
153
154 ! Set default arguments if not present
155 r_min_a = 0.0_dp; if (present(r_min)) r_min_a = r_min
156 coord_a = af_xyz; if (present(coord)) coord_a = coord
157 gb_limit = 4; if (present(mem_limit_gb)) gb_limit = mem_limit_gb
158
159 if (tree%ready) stop "af_init: tree was already initialized"
160 if (n_cell < 2) stop "af_init: n_cell should be >= 2"
161 if (btest(n_cell, 0)) stop "af_init: n_cell should be even"
162 if (gb_limit <= 0) stop "af_init: mem_limit_gb should be > 0"
163 if (coord_a == af_cyl .and. ndim /= 2) stop "af_init: cyl. coords only in 2d"
164 if (tree%n_var_cell <= 0) stop "af_init: no cell-centered variables present"
165
166 do lvl = af_min_lvl, af_max_lvl
167 allocate(tree%lvls(lvl)%ids(0))
168 allocate(tree%lvls(lvl)%leaves(0))
169 allocate(tree%lvls(lvl)%parents(0))
170 end do
171
172 tree%n_cell = n_cell
173 tree%r_base = r_min_a
174 tree%dr_base = (r_max - r_min_a) / grid_size
175 tree%highest_id = 0
176 tree%highest_lvl = 0
177 tree%coord_t = coord_a
178
179 if (present(box_limit)) then
180 if (box_limit <= 0) stop "af_init: box_limit should be > 0"
181 tree%box_limit = box_limit
182 else
183 ! Calculate size of a box
184 box_bytes = af_box_bytes(n_cell, tree%n_var_cell, tree%n_var_face)
185 tree%box_limit = nint(gb_limit * 2.0_dp**30 / box_bytes)
186 end if
187
188 ! Allocate the full list of boxes
189 allocate(tree%boxes(tree%box_limit))
190
191 ! This list can probably be a bit smaller
192 tree%n_removed_ids = 0
193 allocate(tree%removed_ids(tree%box_limit))
194
195 ! Initialize list of cell-centered variables with methods
196 if (.not. allocated(tree%cc_auto_vars)) &
197 allocate(tree%cc_auto_vars(0))
198 if (.not. allocated(tree%cc_func_vars)) &
199 allocate(tree%cc_func_vars(0))
200
201 call af_set_coarse_grid(tree, grid_size, periodic)
202
203 end subroutine af_init
204
205 !> Create the coarse grid
206 subroutine af_set_coarse_grid(tree, coarse_grid_size, periodic_dims)
207 !> Tree for which we set the base
208 type(af_t), intent(inout) :: tree
209 !> Size of coarse grid (in cells)
210 integer, intent(in) :: coarse_grid_size(NDIM)
211 !> Whether dimensions are periodic (default: false)
212 logical, intent(in), optional :: periodic_dims(NDIM)
213 logical :: periodic(NDIM)
214 integer :: nx(NDIM), ix(NDIM), IJK, id, n_boxes, nb
215 integer :: n, iv
216 integer, allocatable :: id_array(DTIMES(:))
217
218 if (tree%highest_id > 0) &
219 error stop "af_set_coarse_grid: this tree already has boxes"
220 if (.not. allocated(tree%boxes)) &
221 error stop "af_set_coarse_grid: tree not initialized"
222 if (any(coarse_grid_size < tree%n_cell)) &
223 error stop "af_set_coarse_grid: coarse_grid_size < tree%n_cell"
224 if (any(modulo(coarse_grid_size, tree%n_cell) /= 0)) &
225 error stop "af_set_coarse_grid: coarse_grid_size not divisible by tree%n_cell"
226
227 periodic(:) = .false.; if (present(periodic_dims)) periodic = periodic_dims
228
229 tree%coarse_grid_size(1:ndim) = coarse_grid_size
230 tree%periodic(1:ndim) = periodic
231 tree%ready = .true.
232
233 nx = coarse_grid_size / tree%n_cell
234 n_boxes = product(nx)
235
236 ! For easy lookup of box neighbors
237 call create_index_array(nx, periodic, id_array)
238
239 ! Check if we have enough space
240 if (n_boxes > size(tree%boxes(:))) &
241 error stop "Not enough memory available for coarse grid"
242
243 ! Create level 1
244 deallocate(tree%lvls(1)%ids)
245 allocate(tree%lvls(1)%ids(n_boxes))
246
247 ! The ids are simply 1, 2, 3, ..., N
248 call get_free_ids(tree, tree%lvls(1)%ids)
249 tree%lvls(1)%leaves = tree%lvls(1)%ids
250
251 ! Loop over the boxes and set their neighbors
252#if NDIM == 1
253 do i = 1, nx(1)
254 id = id_array(ijk)
255 tree%boxes(id)%lvl = 1
256 tree%boxes(id)%ix = [ijk]
257 tree%boxes(id)%dr = tree%dr_base
258 tree%boxes(id)%r_min = tree%r_base + &
259 (tree%boxes(id)%ix - 1) * tree%dr_base * tree%n_cell
260 tree%boxes(id)%n_cell = tree%n_cell
261 tree%boxes(id)%coord_t = tree%coord_t
262
263 tree%boxes(id)%parent = af_no_box
264 tree%boxes(id)%children(:) = af_no_box
265
266 ! Connectivity
267 do nb = 1, af_num_neighbors
268 ix = [ijk] + af_neighb_dix(:, nb)
269 tree%boxes(id)%neighbors(nb) = id_array(ix(1))
270 end do
271 tree%boxes(id)%neighbor_mat = id_array(i-1:i+1)
272
273 call af_init_box(tree, id)
274 end do
275#elif NDIM == 2
276 do j = 1, nx(2)
277 do i = 1, nx(1)
278 id = id_array(ijk)
279 tree%boxes(id)%lvl = 1
280 tree%boxes(id)%ix = [ijk]
281 tree%boxes(id)%dr = tree%dr_base
282 tree%boxes(id)%r_min = tree%r_base + &
283 (tree%boxes(id)%ix - 1) * tree%dr_base * tree%n_cell
284 tree%boxes(id)%n_cell = tree%n_cell
285 tree%boxes(id)%coord_t = tree%coord_t
286
287 tree%boxes(id)%parent = af_no_box
288 tree%boxes(id)%children(:) = af_no_box
289
290 ! Connectivity
291 do nb = 1, af_num_neighbors
292 ix = [ijk] + af_neighb_dix(:, nb)
293 tree%boxes(id)%neighbors(nb) = id_array(ix(1), ix(2))
294 end do
295 tree%boxes(id)%neighbor_mat = id_array(i-1:i+1, j-1:j+1)
296
297 call af_init_box(tree, id)
298 end do
299 end do
300#elif NDIM == 3
301 do k = 1, nx(3)
302 do j = 1, nx(2)
303 do i = 1, nx(1)
304 id = id_array(ijk)
305 tree%boxes(id)%lvl = 1
306 tree%boxes(id)%ix = [ijk]
307 tree%boxes(id)%dr = tree%dr_base
308 tree%boxes(id)%r_min = tree%r_base + &
309 (tree%boxes(id)%ix - 1) * tree%dr_base * tree%n_cell
310 tree%boxes(id)%n_cell = tree%n_cell
311 tree%boxes(id)%coord_t = tree%coord_t
312
313 tree%boxes(id)%parent = af_no_box
314 tree%boxes(id)%children(:) = af_no_box
315
316 ! Connectivity
317 do nb = 1, af_num_neighbors
318 ix = [ijk] + af_neighb_dix(:, nb)
319 tree%boxes(id)%neighbors(nb) = id_array(ix(1), ix(2), ix(3))
320 end do
321 tree%boxes(id)%neighbor_mat = id_array(i-1:i+1, j-1:j+1, k-1:k+1)
322
323 call af_init_box(tree, id)
324 end do
325 end do
326 end do
327#endif
328
329 tree%highest_lvl = 1
330
331 ! Set values for variables with a 'funcval'
332 do i = 1, size(tree%lvls(1)%ids)
333 id = tree%lvls(1)%ids(i)
334 do n = 1, size(tree%cc_func_vars)
335 iv = tree%cc_func_vars(n)
336 call tree%cc_methods(iv)%funcval(tree%boxes(id), iv)
337 end do
338 end do
339
340 end subroutine af_set_coarse_grid
341
342 !> Set the methods for a cell-centered variable
343 subroutine af_set_cc_methods(tree, iv, bc, rb, prolong, restrict, &
344 bc_custom, funcval, prolong_limiter)
345 use m_af_ghostcell, only: af_gc_interp
346 use m_af_prolong, only: af_prolong_linear
347 use m_af_restrict, only: af_restrict_box
348 use m_af_limiters
349 type(af_t), intent(inout) :: tree !< Tree to operate on
350 integer, intent(in) :: iv !< Index of variable
351 procedure(af_subr_bc), optional :: bc !< Boundary condition method
352 procedure(af_subr_rb), optional :: rb !< Refinement boundary method
353 procedure(af_subr_prolong), optional :: prolong !< Prolongation method
354 procedure(af_subr_restrict), optional :: restrict !< Restriction method
355 procedure(af_subr_bc_custom), optional :: bc_custom !< Custom b.c. method
356 procedure(af_subr_funcval), optional :: funcval !< Variable defined by function
357 !< Type of limiter to use for prolongation (of values or ghost cells)
358 integer, intent(in), optional :: prolong_limiter
359 integer :: i
360
361 if (tree%has_cc_method(iv)) then
362 print *, "Cannot call af_set_cc_methods twice for ", &
363 trim(tree%cc_names(iv))
364 error stop
365 end if
366
367 ! Set methods for the variable and its copies
368 do i = iv, iv + tree%cc_num_copies(iv) - 1
369 if (present(bc)) then
370 tree%cc_methods(i)%bc => bc
371 else if (present(bc_custom)) then
372 tree%cc_methods(i)%bc_custom => bc_custom
373 else if (.not. present(funcval)) then
374 error stop "af_set_cc_methods: bc, bc_custom or funcval required"
375 end if
376
377 if (present(funcval)) then
378 tree%cc_methods(i)%funcval => funcval
379 end if
380
381 if (present(rb)) then
382 tree%cc_methods(i)%rb => rb
383 else
384 tree%cc_methods(i)%rb => af_gc_interp
385 end if
386
387 if (present(prolong)) then
388 tree%cc_methods(i)%prolong => prolong
389 else
390 tree%cc_methods(i)%prolong => af_prolong_linear
391 end if
392
393 if (present(restrict)) then
394 tree%cc_methods(i)%restrict => restrict
395 else
396 tree%cc_methods(i)%restrict => af_restrict_box
397 end if
398
399 if (present(prolong_limiter)) then
400 tree%cc_methods(i)%prolong_limiter = prolong_limiter
401 else if (ndim < 3) then
402 tree%cc_methods(i)%prolong_limiter = af_limiter_mc_t
403 else
404 ! To ensure the interpolation of ghost cells near refinement
405 ! boundaries is non-negative and does not create new maxima in 3D,
406 ! this limiter can be used
407 tree%cc_methods(i)%prolong_limiter = af_limiter_gminmod43_t
408 end if
409
410 tree%has_cc_method(i) = .true.
411 end do
412
413 if (.not. allocated(tree%cc_auto_vars)) &
414 allocate(tree%cc_auto_vars(0))
415 if (.not. allocated(tree%cc_func_vars)) &
416 allocate(tree%cc_func_vars(0))
417
418
419 ! Append only original variable, so that the copies are not automatically
420 ! prolongated etc.
421 if (present(funcval)) then
422 tree%cc_func_vars = [tree%cc_func_vars, iv]
423 else
424 tree%cc_auto_vars = [tree%cc_auto_vars, iv]
425 end if
426
427 end subroutine af_set_cc_methods
428
429 !> "Destroy" the data in a tree. Since we don't use pointers, you can also
430 !> just let a tree get out of scope
431 subroutine af_destroy(tree)
432 type(af_t), intent(out) :: tree
433 end subroutine af_destroy
434
435 !> Create an array for easy lookup of indices
436 subroutine create_index_array(nx, periodic, id_array)
437 integer, intent(in) :: nx(NDIM)
438 logical, intent(in) :: periodic(NDIM)
439 integer, intent(inout), allocatable :: id_array(DTIMES(:))
440 integer :: IJK
441
442#if NDIM == 1
443 allocate(id_array(0:nx(1)+1))
444#elif NDIM == 2
445 allocate(id_array(0:nx(1)+1, 0:nx(2)+1))
446#elif NDIM == 3
447 allocate(id_array(0:nx(1)+1, 0:nx(2)+1, 0:nx(3)+1))
448#endif
449
450 id_array = af_phys_boundary
451
452#if NDIM == 1
453 do i = 1, nx(1)
454 id_array(i) = i
455 end do
456
457 if (periodic(1)) then
458 id_array(0) = id_array(nx(1))
459 id_array(nx(1)+1) = id_array(1)
460 end if
461#elif NDIM == 2
462 do j = 1, nx(2)
463 do i = 1, nx(1)
464 id_array(i, j) = (j-1) * nx(1) + i
465 end do
466 end do
467
468 if (periodic(1)) then
469 id_array(0, :) = id_array(nx(1), :)
470 id_array(nx(1)+1, :) = id_array(1, :)
471 end if
472
473 if (periodic(2)) then
474 id_array(:, 0) = id_array(:, nx(2))
475 id_array(:, nx(2)+1) = id_array(:, 1)
476 end if
477#elif NDIM == 3
478 do k = 1, nx(3)
479 do j = 1, nx(2)
480 do i = 1, nx(1)
481 id_array(i, j, k) = (k-1) * nx(2) * nx(1) + (j-1) * nx(1) + i
482 end do
483 end do
484 end do
485
486 if (periodic(1)) then
487 id_array(0, :, :) = id_array(nx(1), :, :)
488 id_array(nx(1)+1, :, :) = id_array(1, :, :)
489 end if
490
491 if (periodic(2)) then
492 id_array(:, 0, :) = id_array(:, nx(2), :)
493 id_array(:, nx(2)+1, :) = id_array(:, 1, :)
494 end if
495
496 if (periodic(3)) then
497 id_array(:, :, 0) = id_array(:, :, nx(3))
498 id_array(:, :, nx(3)+1) = id_array(:, :, 1)
499 end if
500#endif
501 end subroutine create_index_array
502
503 !> Create a list of leaves and a list of parents for a level
504 subroutine set_leaves_parents(boxes, level)
505 type(box_t), intent(in) :: boxes(:) !< List of boxes
506 type(lvl_t), intent(inout) :: level !< Level type which contains the indices of boxes
507 integer :: i, id, i_leaf, i_parent
508 integer :: n_parents, n_leaves
509
510 n_parents = count(af_has_children(boxes(level%ids)))
511 n_leaves = size(level%ids) - n_parents
512
513 if (n_parents /= size(level%parents)) then
514 deallocate(level%parents)
515 allocate(level%parents(n_parents))
516 end if
517
518 if (n_leaves /= size(level%leaves)) then
519 deallocate(level%leaves)
520 allocate(level%leaves(n_leaves))
521 end if
522
523 i_leaf = 0
524 i_parent = 0
525 do i = 1, size(level%ids)
526 id = level%ids(i)
527 if (af_has_children(boxes(id))) then
528 i_parent = i_parent + 1
529 level%parents(i_parent) = id
530 else
531 i_leaf = i_leaf + 1
532 level%leaves(i_leaf) = id
533 end if
534 end do
535 end subroutine set_leaves_parents
536
537 !> Mark box as active and allocate data storage for a box, for its cell- and
538 !> face-centered data
539 subroutine af_init_box(tree, id)
540 type(af_t), intent(inout) :: tree !< Tree
541 integer, intent(in) :: id !< Box id
542 integer :: nc, ix, nb
543 logical :: new_box
544
545 associate(box => tree%boxes(id))
546 nc = tree%n_cell
547 box%in_use = .true.
548 new_box = .not. allocated(box%cc)
549
550 if (new_box) then
551 allocate(box%cc(dtimes(0:nc+1), tree%n_var_cell))
552 allocate(box%fc(dtimes(nc+1), ndim, tree%n_var_face))
553 end if
554
555 ! Initialize to zero
556 box%cc = 0
557 box%fc = 0
558
559 ! Allocate storage for boundary conditions
560 box%n_bc = count(box%neighbors < af_no_box)
561 allocate(box%bc_index_to_nb(box%n_bc))
562 allocate(box%bc_coords(ndim, nc**(ndim-1), box%n_bc))
563 allocate(box%bc_val(nc**(ndim-1), tree%n_var_cell, box%n_bc))
564 allocate(box%bc_type(tree%n_var_cell, box%n_bc))
565 box%bc_val = 0
566 box%bc_type = 0
567
568 ! Set face coordinates
569 ix = 0
570 do nb = 1, af_num_neighbors
571 if (box%neighbors(nb) < af_no_box) then
572 ix = ix + 1
573 box%bc_index_to_nb(ix) = nb
574 box%nb_to_bc_index(nb) = ix
575 call af_get_face_coords(box, nb, box%bc_coords(:, :, ix))
576 end if
577 end do
578 end associate
579 end subroutine af_init_box
580
581 !> Mark box as inactive, but keep storage for cell- and face-centered data to
582 !> avoid reallocating this
583 subroutine af_deactivate_box(box)
584 type(box_t), intent(inout) :: box
585
586 box%in_use = .false.
587 box%tag = af_init_tag
588 box%n_stencils = 0
589 if (allocated(box%stencils)) deallocate(box%stencils)
590 deallocate(box%bc_index_to_nb, box%bc_coords, &
591 box%bc_val, box%bc_type)
592 end subroutine af_deactivate_box
593
594 ! Set the neighbors of id (using their parent)
595 subroutine set_neighbs(boxes, id)
596 type(box_t), intent(inout) :: boxes(:)
597 integer, intent(in) :: id
598 integer :: nb, nb_id, IJK
599
600 do kji_do(-1, 1)
601 if (boxes(id)%neighbor_mat(ijk) == af_no_box) then
602 nb_id = find_neighb(boxes, id, [ijk])
603 if (nb_id > af_no_box) then
604 boxes(id)%neighbor_mat(ijk) = nb_id
605#if NDIM == 1
606 boxes(nb_id)%neighbor_mat(-i) = id
607#elif NDIM == 2
608 boxes(nb_id)%neighbor_mat(-i, -j) = id
609#elif NDIM == 3
610 boxes(nb_id)%neighbor_mat(-i, -j, -k) = id
611#endif
612 end if
613 end if
614 end do; close_do
615
616 do nb = 1, af_num_neighbors
617 if (boxes(id)%neighbors(nb) == af_no_box) then
618#if NDIM == 1
619 nb_id = boxes(id)%neighbor_mat(af_neighb_dix(1, nb))
620#elif NDIM == 2
621 nb_id = boxes(id)%neighbor_mat(af_neighb_dix(1, nb), &
622 af_neighb_dix(2, nb))
623#elif NDIM == 3
624 nb_id = boxes(id)%neighbor_mat(af_neighb_dix(1, nb), &
625 af_neighb_dix(2, nb), af_neighb_dix(3, nb))
626#endif
627 if (nb_id > af_no_box) then
628 boxes(id)%neighbors(nb) = nb_id
629 boxes(nb_id)%neighbors(af_neighb_rev(nb)) = id
630 end if
631 end if
632 end do
633 end subroutine set_neighbs
634
635 !> Get the id of all neighbors of boxes(id), through its parent
636 function find_neighb(boxes, id, dix) result(nb_id)
637 type(box_t), intent(in) :: boxes(:) !< List with all the boxes
638 integer, intent(in) :: id !< Box whose neighbor we are looking for
639 integer, intent(in) :: dix(ndim)
640 integer :: nb_id, p_id, c_ix, dix_c(ndim)
641
642 p_id = boxes(id)%parent
643 c_ix = af_ix_to_ichild(boxes(id)%ix)
644
645 ! Check if neighbor is in same direction as dix is (low/high). If so, use
646 ! neighbor of parent
647 where ((dix == -1) .eqv. af_child_low(:, c_ix))
648 dix_c = dix
649 elsewhere
650 dix_c = 0
651 end where
652
653 p_id = boxes(p_id)%neighbor_mat(dindex(dix_c))
654
655 if (p_id <= af_no_box) then
656 nb_id = p_id
657 else
658 c_ix = af_ix_to_ichild(boxes(id)%ix + dix)
659 nb_id = boxes(p_id)%children(c_ix)
660 end if
661 end function find_neighb
662
663 !> Refine a tree up to a given refinement lvl
664 subroutine af_refine_up_to_lvl(tree, lvl)
665 type(af_t), intent(inout) :: tree !< The tree to adjust
666 integer, intent(in) :: lvl !< Refine up to this lvl
667 type(ref_info_t) :: ref_info
668
669 if (lvl < tree%highest_lvl) error stop "tree already above level"
670
671 do
672 call af_adjust_refinement(tree, ref_routine, ref_info)
673 if (ref_info%n_add == 0) exit
674 end do
675
676 contains
677
678 subroutine ref_routine(box, cell_flags)
679 type(box_t), intent(in) :: box
680 integer, intent(out) :: cell_flags(dtimes(box%n_cell))
681
682 if (box%lvl < lvl) then
683 cell_flags = af_do_ref
684 else
685 cell_flags = af_keep_ref
686 end if
687 end subroutine ref_routine
688
689 end subroutine af_refine_up_to_lvl
690
691 !> Adjust the refinement of a tree using the user-supplied ref_subr. The
692 !> optional argument ref_buffer controls over how many cells neighbors are
693 !> affected by refinement flags.
694 !>
695 !> On input, the tree should be balanced. On output, the tree is still
696 !> balanced, and its refinement is updated (with at most one level per call).
697 subroutine af_adjust_refinement(tree, ref_subr, ref_info, ref_buffer, &
698 ref_links)
699 type(af_t), intent(inout) :: tree !< The tree to adjust
700 procedure(af_subr_ref) :: ref_subr !< Refinement function
701 type(ref_info_t), intent(inout) :: ref_info !< Information about refinement
702 integer, intent(in), optional :: ref_buffer !< Buffer width (in cells)
703 !> Lists of linked boxes which should have the same refinement
704 integer, intent(in), optional :: ref_links(:, :)
705 integer :: lvl, id, i, c_ids(af_num_children), i_ch
706 integer :: i_add, i_rm, n_ch, n_add
707 integer, allocatable :: ref_flags(:)
708 integer :: ref_buffer_val
709
710 if (.not. tree%ready) stop "Tree not ready"
711
712 ref_buffer_val = 0 ! Default buffer width (in cells) around refinement
713 if (present(ref_buffer)) ref_buffer_val = ref_buffer
714
715 if (ref_buffer_val < 0) &
716 error stop "af_adjust_refinement: ref_buffer < 0"
717 if (ref_buffer_val > tree%n_cell) &
718 error stop "af_adjust_refinement: ref_buffer > tree%n_cell"
719
720 allocate(ref_flags(tree%highest_id))
721
722 ! Set refinement values for all boxes. Only two flags are used below:
723 ! af_refine and af_derefine. Other values are ignored.
724 call consistent_ref_flags(tree, ref_flags, ref_subr, &
725 ref_buffer_val, ref_links)
726
727 ! To store ids of removed boxes
728 n_ch = af_num_children
729 ref_info%n_rm = n_ch * count(ref_flags == af_derefine)
730 if (allocated(ref_info%rm)) deallocate(ref_info%rm)
731 allocate(ref_info%rm(ref_info%n_rm))
732
733 ! To store ids of new boxes per level
734 ref_info%n_add = n_ch * count(ref_flags == af_refine)
735 if (allocated(ref_info%lvls)) deallocate(ref_info%lvls)
736 allocate(ref_info%lvls(tree%highest_lvl+1))
737
738 ! There can be no new children at level 1
739 allocate(ref_info%lvls(1)%add(0))
740
741 do lvl = 1, tree%highest_lvl
742 ! Number of newly added boxes to the next level
743 n_add = n_ch * count(ref_flags(tree%lvls(lvl)%ids) == af_refine)
744 allocate(ref_info%lvls(lvl+1)%add(n_add))
745 end do
746
747 i_rm = 0
748
749 do lvl = 1, af_max_lvl-1 ! The loop exits when it encounters an empty level
750 i_add = 0
751
752 do i = 1, size(tree%lvls(lvl)%ids)
753 id = tree%lvls(lvl)%ids(i)
754
755 if (id > size(ref_flags)) then
756 cycle ! This is a newly added box
757 else if (ref_flags(id) == af_refine) then
758 ! Add children. First need to get num_children free id's
759 call get_free_ids(tree, c_ids)
760 call add_children(tree, id, c_ids)
761 ref_info%lvls(lvl+1)%add(i_add+1:i_add+n_ch) = &
762 tree%boxes(id)%children
763 i_add = i_add + n_ch
764 else if (ref_flags(id) == af_derefine) then
765 ! Remove children
766 call auto_restrict(tree, id)
767 ref_info%rm(i_rm+1:i_rm+n_ch) = tree%boxes(id)%children
768 i_rm = i_rm + n_ch
769 call remove_children(tree, id)
770 end if
771 end do
772
773 ! Update leaves / parents
774 call set_leaves_parents(tree%boxes, tree%lvls(lvl))
775
776 ! Set next level ids to children of this level
777 call set_child_ids(tree%lvls(lvl)%parents, &
778 tree%lvls(lvl+1)%ids, tree%boxes)
779
780 ! Update connectivity of new children
781 do i = 1, size(tree%lvls(lvl)%parents)
782 id = tree%lvls(lvl)%parents(i)
783 if (ref_flags(id) == af_refine) then
784 do i_ch = 1, af_num_children
785 call set_neighbs(tree%boxes, tree%boxes(id)%children(i_ch))
786 end do
787 end if
788 end do
789
790 if (size(tree%lvls(lvl+1)%ids) == 0) exit
791 end do
792
793 tree%highest_lvl = lvl
794
795 ! Update the list of removed boxes. We do this at the end so that they are
796 ! not re-used in the loop above.
797 i = tree%n_removed_ids
798 tree%removed_ids(i+1:i+ref_info%n_rm) = ref_info%rm(:)
799 tree%n_removed_ids = tree%n_removed_ids + ref_info%n_rm
800
801 ! Set the highest id in use
802 do id = tree%highest_id, 1, -1
803 if (tree%boxes(id)%in_use) exit
804 end do
805 tree%highest_id = id
806
807 ! Clean up list of removed boxes
808 i_rm = tree%n_removed_ids
809 i = count(tree%removed_ids(1:i_rm) <= tree%highest_id)
810 tree%removed_ids(1:i) = pack(tree%removed_ids(1:i_rm), &
811 mask=tree%removed_ids(1:i_rm) <= tree%highest_id)
812 tree%n_removed_ids = i
813
814 ! We still have to update leaves and parents for the last level, which is
815 ! either lvl+1 or af_max_lvl. Note that lvl+1 is empty now, but maybe it was
816 ! not not empty before, and that af_max_lvl is skipped in the above loop.
817 lvl = min(lvl+1, af_max_lvl)
818 call set_leaves_parents(tree%boxes, tree%lvls(lvl))
819
820 call auto_prolong(tree, ref_info)
821
822 end subroutine af_adjust_refinement
823
824 !> Try to automatically restrict to box with index id
825 subroutine auto_restrict(tree, id)
826 type(af_t), intent(inout) :: tree
827 integer, intent(in) :: id
828 integer :: i, iv, i_ch, ch_id
829
830 if (.not. any(tree%has_cc_method(:))) return
831
832 do i_ch = 1, af_num_children
833 ch_id = tree%boxes(id)%children(i_ch)
834 do i = 1, size(tree%cc_auto_vars)
835 iv = tree%cc_auto_vars(i)
836 call tree%cc_methods(iv)%restrict(tree%boxes(ch_id), &
837 tree%boxes(id), [iv])
838 end do
839 end do
840 end subroutine auto_restrict
841
842 !> Try to automatically prolong to all new boxes
843 subroutine auto_prolong(tree, ref_info)
844 use m_af_ghostcell, only: af_gc_box
845 type(af_t), intent(inout) :: tree
846 type(ref_info_t), intent(in) :: ref_info
847 integer :: lvl, i, n, iv, id, p_id
848
849 ! Skip this routine when it won't do anything
850 if (.not. any(tree%has_cc_method(:)) .or. ref_info%n_add == 0) then
851 return
852 end if
853
854 !$omp parallel private(lvl, i, n, iv, id, p_id)
855 do lvl = 1, tree%highest_lvl
856 !$omp do
857 do i = 1, size(ref_info%lvls(lvl)%add)
858 id = ref_info%lvls(lvl)%add(i)
859 p_id = tree%boxes(id)%parent
860
861 do n = 1, size(tree%cc_auto_vars)
862 iv = tree%cc_auto_vars(n)
863 call tree%cc_methods(iv)%prolong(tree%boxes(p_id), &
864 tree%boxes(id), iv, limiter=tree%cc_methods(iv)%prolong_limiter)
865 end do
866 do n = 1, size(tree%cc_func_vars)
867 iv = tree%cc_func_vars(n)
868 call tree%cc_methods(iv)%funcval(tree%boxes(id), iv)
869 end do
870 end do
871 !$omp end do
872
873 !$omp do
874 do i = 1, size(ref_info%lvls(lvl)%add)
875 id = ref_info%lvls(lvl)%add(i)
876 call af_gc_box(tree, id, [tree%cc_auto_vars])
877 end do
878 !$omp end do
879 end do
880 !$omp end parallel
881 end subroutine auto_prolong
882
883 !> Get free ids from the boxes(:) array to store new boxes in. These ids are
884 !> always consecutive.
885 subroutine get_free_ids(tree, ids)
886 type(af_t), intent(inout) :: tree
887 integer, intent(out) :: ids(:) !< Array which will be filled with free box ids
888 integer :: i, highest_id_prev, n_ids
889
890 n_ids = size(ids)
891
892 !> @todo when doing AMR in parallel, perhaps move some of the code outside
893 !> the critical construct
894
895 !$omp critical (crit_free_ids)
896 if (n_ids <= tree%n_removed_ids) then
897 ! Re-use removed boxes
898 do i = 1, n_ids
899 ids(i) = tree%removed_ids(tree%n_removed_ids-n_ids+i)
900 end do
901 tree%n_removed_ids = tree%n_removed_ids - n_ids
902 else
903 ! Add new boxes at the end of the list
904 highest_id_prev = tree%highest_id
905 tree%highest_id = tree%highest_id + n_ids
906
907 if (tree%highest_id > size(tree%boxes)) then
908 print *, "get_free_ids: exceeding memory limit"
909 write(*, '(A,E12.2)') " memory_limit (GByte): ", &
910 tree%box_limit * 0.5_dp**30 * &
911 af_box_bytes(tree%n_cell, tree%n_var_cell, tree%n_var_face)
912 print *, "memory_limit (boxes): ", tree%box_limit
913 print *, "You can increase the memory limit in your call to af_init"
914 print *, "by setting mem_limit_gb to a higher value (in GBytes)"
915 error stop
916 end if
917
918 ids = [(highest_id_prev + i, i=1,n_ids)]
919 end if
920 !$omp end critical (crit_free_ids)
921
922 end subroutine get_free_ids
923
924 !> Given the refinement function, return consistent refinement flags, that
925 !> ensure that the tree is still balanced. Furthermore, it cannot derefine the
926 !> base level, and it cannot refine above af_max_lvl. The argument
927 !> ref_flags is changed: for boxes that will be refined it holds af_refine,
928 !> for boxes that will be derefined it holds af_derefine
929 subroutine consistent_ref_flags(tree, ref_flags, ref_subr, &
930 ref_buffer, ref_links)
931 use omp_lib, only: omp_get_max_threads, omp_get_thread_num
932 type(af_t), intent(inout) :: tree !< Tree for which we set refinement flags
933 integer, intent(inout) :: ref_flags(:) !< List of refinement flags for all boxes(:)
934 procedure(af_subr_ref) :: ref_subr !< User-supplied refinement function.
935 integer, intent(in) :: ref_buffer !< Buffer width (in cells)
936 !> Lists of linked boxes which should have the same refinement
937 integer, intent(in), optional :: ref_links(:, :)
938 integer :: lvl, i, i_ch, ch_id, id
939 integer :: p_id
940 integer :: thread_id
941 integer, allocatable :: tmp_flags(:, :)
942 integer :: cell_flags(DTIMES(tree%n_cell))
943 integer, parameter :: unset_flag = -huge(1)
944
945 ! Set refinement flags for each thread individually, because we sometimes
946 ! modify the refinement flags of neighbors
947 allocate(tmp_flags(size(ref_flags), omp_get_max_threads()))
948
949 tmp_flags(:, :) = unset_flag
950
951 ! Set refinement flags on all leaves and their immediate parents (on other
952 ! boxes the flags would not matter)
953
954 !$omp parallel private(lvl, i, id, p_id, cell_flags, thread_id, i_ch, ch_id)
955 thread_id = omp_get_thread_num() + 1
956
957 do lvl = 1, tree%highest_lvl
958 !$omp do
959 do i = 1, size(tree%lvls(lvl)%leaves)
960 id = tree%lvls(lvl)%leaves(i)
961
962 call ref_subr(tree%boxes(id), cell_flags)
963 call cell_to_ref_flags(cell_flags, tree%n_cell, &
964 tmp_flags(:, thread_id), tree, id, ref_buffer)
965
966 ! If the parent exists, and this is the first child which is itself
967 ! not refined, set refinement flags for the parent
968 if (tree%boxes(id)%lvl > 1) then
969 p_id = tree%boxes(id)%parent
970 do i_ch = 1, af_ix_to_ichild(tree%boxes(id)%ix)-1
971 ch_id = tree%boxes(p_id)%children(i_ch)
972 if (.not. af_has_children(tree%boxes(ch_id))) exit
973 end do
974
975 if (i_ch == af_ix_to_ichild(tree%boxes(id)%ix)) then
976 ! This is the first child which is itself not refined
977 call ref_subr(tree%boxes(p_id), cell_flags)
978 call cell_to_ref_flags(cell_flags, tree%n_cell, &
979 tmp_flags(:, thread_id), tree, p_id, ref_buffer)
980 end if
981 end if
982 end do
983 !$omp end do
984 end do
985 !$omp end parallel
986
987 ! Take the highest value over the threads
988 do i = 1, size(ref_flags)
989 ref_flags(i) = maxval(tmp_flags(i, :))
990 if (ref_flags(i) == unset_flag) ref_flags(i) = af_keep_ref
991 end do
992
993 if (maxval(ref_flags) > af_do_ref .or. minval(ref_flags) < af_rm_ref) &
994 stop "af_adjust_refinement: invalid refinement flag given"
995
996 ! Cannot refine beyond max level
997 do i = 1, size(tree%lvls(af_max_lvl)%ids)
998 id = tree%lvls(af_max_lvl)%ids(i)
999 if (ref_flags(id) == af_do_ref) ref_flags(id) = af_keep_ref
1000 end do
1001
1002 call ensure_two_one_balance(tree, ref_flags)
1003 call handle_derefinement_flags(tree, ref_flags)
1004
1005 if (present(ref_links)) then
1006 do i = 1, size(ref_links, 2)
1007 ref_flags(ref_links(:, i)) = maxval(ref_flags(ref_links(:, i)))
1008 end do
1009 call ensure_two_one_balance(tree, ref_flags)
1010 call handle_derefinement_flags(tree, ref_flags)
1011 end if
1012
1013 end subroutine consistent_ref_flags
1014
1015 !> Adjust refinement flags to ensure 2-1 balance is maintained
1016 subroutine ensure_two_one_balance(tree, ref_flags)
1017 type(af_t), intent(inout) :: tree
1018 integer, intent(inout) :: ref_flags(:)
1019 integer :: lvl, i, id, nb, nb_id
1020 integer :: p_id, p_nb_id
1021
1022 ! Ensure 2-1 balance
1023 do lvl = tree%highest_lvl, 1, -1
1024 do i = 1, size(tree%lvls(lvl)%leaves) ! We only check leaf tree%boxes
1025 id = tree%lvls(lvl)%leaves(i)
1026
1027 if (ref_flags(id) == af_do_ref .or. ref_flags(id) == af_refine) then
1028 ref_flags(id) = af_refine ! Mark for actual refinement
1029
1030 ! Ensure we will have the necessary neighbors
1031 do nb = 1, af_num_neighbors
1032 nb_id = tree%boxes(id)%neighbors(nb)
1033 if (nb_id == af_no_box) then
1034 ! Mark the parent containing neighbor for refinement
1035 p_id = tree%boxes(id)%parent
1036 p_nb_id = tree%boxes(p_id)%neighbors(nb)
1037 ref_flags(p_nb_id) = af_refine
1038 end if
1039 end do
1040
1041 else if (ref_flags(id) == af_rm_ref) then
1042 ! Ensure we do not remove a required neighbor
1043 do nb = 1, af_num_neighbors
1044 nb_id = tree%boxes(id)%neighbors(nb)
1045 if (nb_id > af_no_box) then
1046 if (af_has_children(tree%boxes(nb_id)) .or. &
1047 ref_flags(nb_id) > af_keep_ref) then
1048 ref_flags(id) = af_keep_ref
1049 exit
1050 end if
1051 end if
1052 end do
1053 end if
1054
1055 end do
1056 end do
1057 end subroutine ensure_two_one_balance
1058
1059 subroutine handle_derefinement_flags(tree, ref_flags)
1060 type(af_t), intent(inout) :: tree
1061 integer, intent(inout) :: ref_flags(:)
1062 integer :: lvl, i, id, c_ids(af_num_children)
1063
1064 ! Make the (de)refinement flags consistent for blocks with children
1065 do lvl = tree%highest_lvl-1, 1, -1
1066 do i = 1, size(tree%lvls(lvl)%parents)
1067 id = tree%lvls(lvl)%parents(i)
1068 c_ids = tree%boxes(id)%children
1069
1070 ! Only consider boxes for which at least one child is a leaf
1071 if (all(af_has_children(tree%boxes(c_ids)))) cycle
1072
1073 ! Can only remove children if they are all marked for
1074 ! derefinement, and the box itself not for refinement.
1075 if (all(ref_flags(c_ids) == af_rm_ref) .and. &
1076 ref_flags(id) <= af_keep_ref) then
1077 ref_flags(id) = af_derefine
1078 else
1079 ref_flags(id) = af_keep_ref
1080 ! The children cannot be removed. This information is useful when
1081 ! the modify_refinement() routine is used. Make sure not to
1082 ! override previously set derefinement flags.
1083 where (ref_flags(c_ids) /= af_derefine)
1084 ref_flags(c_ids) = max(ref_flags(c_ids), af_keep_ref)
1085 end where
1086 end if
1087 end do
1088 end do
1089
1090 end subroutine handle_derefinement_flags
1091
1092 !> Given the cell refinement flags of a box, set the refinement flag for that
1093 !> box and potentially also its neighbors (in case of refinement near a
1094 !> boundary).
1095 subroutine cell_to_ref_flags(cell_flags, nc, ref_flags, tree, id, &
1096 ref_buffer)
1097 use m_af_utils, only: af_get_loc
1098 integer, intent(in) :: nc !< n_cell for the box
1099 integer, intent(in) :: cell_flags(DTIMES(nc)) !< Cell refinement flags
1100 integer, intent(inout) :: ref_flags(:) !< Box refinement flags for this thread
1101 type(af_t), intent(in) :: tree !< Full tree
1102 integer, intent(in) :: id !< Which box is considered
1103 integer, intent(in) :: ref_buffer !< Buffer cells around refinement
1104 integer :: ix0(NDIM), ix1(NDIM), IJK, nb_id
1105
1106 if (minval(cell_flags) < af_rm_ref .or. &
1107 maxval(cell_flags) > af_do_ref) then
1108 error stop "Error: invalid cell flags given"
1109 end if
1110
1111 ! Check whether the box needs to be refined or keep its refinement
1112 if (any(cell_flags == af_do_ref)) then
1113 ref_flags(id) = af_do_ref
1114 else if (any(cell_flags == af_keep_ref)) then
1115 ref_flags(id) = max(ref_flags(id), af_keep_ref)
1116 else ! All flags are af_rm_ref
1117 ref_flags(id) = max(ref_flags(id), af_rm_ref)
1118 end if
1119
1120 if (ref_buffer <= 0) return ! No need to check neighbors
1121
1122 ! Check whether neighbors also require refinement, which happens when cells
1123 ! close to the neighbor are flagged.
1124 do kji_do(-1,1)
1125 if (all([ijk] == 0)) cycle
1126
1127 nb_id = tree%boxes(id)%neighbor_mat(ijk)
1128
1129 ! Skip neighbors that are not there
1130 if (nb_id <= af_no_box) cycle
1131
1132 ! Compute index range relevant for neighbor
1133 ix0 = 1
1134 ix1 = nc
1135 where ([ijk] == 1)
1136 ix0 = nc - ref_buffer + 1
1137 ix1 = nc
1138 elsewhere ([ijk] == -1)
1139 ix0 = 1
1140 ix1 = ref_buffer
1141 end where
1142
1143 if (any(cell_flags(dslice(ix0, ix1)) == af_do_ref)) then
1144 ref_flags(nb_id) = af_do_ref
1145 end if
1146 end do; close_do
1147
1148 end subroutine cell_to_ref_flags
1149
1150 !> Remove the children of box id
1151 subroutine remove_children(tree, id)
1152 type(af_t), intent(inout) :: tree
1153 integer, intent(in) :: id !< Id of box whose children will be removed
1154 integer :: ic, c_id, nb_id, nb_rev, nb, IJK
1155
1156 do ic = 1, af_num_children
1157 c_id = tree%boxes(id)%children(ic)
1158
1159 ! Remove from neighbors
1160 do nb = 1, af_num_neighbors
1161 nb_id = tree%boxes(c_id)%neighbors(nb)
1162 if (nb_id > af_no_box) then
1163 nb_rev = af_neighb_rev(nb)
1164 tree%boxes(nb_id)%neighbors(nb_rev) = af_no_box
1165 end if
1166 end do
1167
1168 do kji_do(-1,1)
1169 nb_id = tree%boxes(c_id)%neighbor_mat(ijk)
1170 if (nb_id > af_no_box) then
1171#if NDIM == 1
1172 tree%boxes(nb_id)%neighbor_mat(-i) = af_no_box
1173#elif NDIM == 2
1174 tree%boxes(nb_id)%neighbor_mat(-i, -j) = af_no_box
1175#elif NDIM == 3
1176 tree%boxes(nb_id)%neighbor_mat(-i, -j, -k) = af_no_box
1177#endif
1178 end if
1179 end do; close_do
1180
1181 call af_deactivate_box(tree%boxes(c_id))
1182 end do
1183
1184 tree%boxes(id)%children = af_no_box
1185 end subroutine remove_children
1186
1187 !> Add children to box id, using the indices in c_ids
1188 subroutine add_children(tree, id, c_ids)
1189 type(af_t), intent(inout) :: tree !< Tree
1190 integer, intent(in) :: id !< Id of box that gets children
1191 integer, intent(in) :: c_ids(af_num_children) !< Free ids for the children
1192 integer :: i, nb, child_nb(2**(NDIM-1))
1193 integer :: c_id, c_ix_base(NDIM), dix(NDIM)
1194
1195 associate(boxes => tree%boxes)
1196 boxes(id)%children = c_ids
1197 c_ix_base = 2 * boxes(id)%ix - 1
1198
1199 do i = 1, af_num_children
1200 c_id = c_ids(i)
1201 boxes(c_id)%ix = c_ix_base + af_child_dix(:,i)
1202 boxes(c_id)%lvl = boxes(id)%lvl+1
1203 boxes(c_id)%parent = id
1204 boxes(c_id)%tag = af_init_tag
1205 boxes(c_id)%children = af_no_box
1206 boxes(c_id)%neighbors = af_no_box
1207 boxes(c_id)%neighbor_mat = af_no_box
1208 boxes(c_id)%neighbor_mat(dtimes(0)) = c_id
1209 boxes(c_id)%n_cell = boxes(id)%n_cell
1210 boxes(c_id)%coord_t = boxes(id)%coord_t
1211 boxes(c_id)%dr = 0.5_dp * boxes(id)%dr
1212 boxes(c_id)%r_min = boxes(id)%r_min + 0.5_dp * boxes(id)%dr * &
1213 af_child_dix(:,i) * boxes(id)%n_cell
1214 end do
1215
1216 ! Set boundary conditions at children
1217 do nb = 1, af_num_neighbors
1218 if (boxes(id)%neighbors(nb) < af_no_box) then
1219 child_nb = c_ids(af_child_adj_nb(:, nb)) ! Neighboring children
1220 boxes(child_nb)%neighbors(nb) = boxes(id)%neighbors(nb)
1221 dix = af_neighb_dix(:, nb)
1222 boxes(child_nb)%neighbor_mat(dindex(dix)) = &
1223 boxes(id)%neighbors(nb)
1224 end if
1225 end do
1226 end associate
1227
1228 ! Have to call this after setting boundary conditions
1229 do i = 1, af_num_children
1230 call af_init_box(tree, c_ids(i))
1231 end do
1232
1233 end subroutine add_children
1234
1235 !> Create a list c_ids(:) of all the children of p_ids(:). This is used after
1236 !> a level has been refined.
1237 subroutine set_child_ids(p_ids, c_ids, boxes)
1238 integer, intent(in) :: p_ids(:) !< All the parents ids
1239 integer, allocatable, intent(inout) :: c_ids(:) !< Output: all the children's ids
1240 type(box_t), intent(in) :: boxes(:) !< List of all the boxes
1241 integer :: i, i0, i1, n_children
1242
1243 n_children = af_num_children * size(p_ids)
1244 if (n_children /= size(c_ids)) then
1245 deallocate(c_ids)
1246 allocate(c_ids(n_children))
1247 end if
1248
1249 do i = 1, size(p_ids)
1250 i1 = i * af_num_children
1251 i0 = i1 - af_num_children + 1
1252 c_ids(i0:i1) = boxes(p_ids(i))%children
1253 end do
1254 end subroutine set_child_ids
1255
1256 !> Restrict fluxes from children to parents on refinement boundaries.
1257 subroutine af_consistent_fluxes(tree, f_ixs)
1258 type(af_t), intent(inout) :: tree !< Tree to operate on
1259 integer, intent(in) :: f_ixs(:) !< Indices of the fluxes
1260 integer :: lvl, i, id, nb, nb_id
1261
1262 if (.not. tree%ready) stop "Tree not ready"
1263 !$omp parallel private(lvl, i, id, nb, nb_id)
1264 do lvl = 1, tree%highest_lvl-1
1265 !$omp do
1266 do i = 1, size(tree%lvls(lvl)%parents)
1267 id = tree%lvls(lvl)%parents(i)
1268 do nb = 1, af_num_neighbors
1269 nb_id = tree%boxes(id)%neighbors(nb)
1270
1271 ! If the neighbor exists and has no children, set flux
1272 if (nb_id > af_no_box) then
1273 if (.not. af_has_children(tree%boxes(nb_id))) then
1274 call flux_from_children(tree%boxes, id, nb, f_ixs)
1275 end if
1276 end if
1277 end do
1278 end do
1279 !$omp end do
1280 end do
1281 !$omp end parallel
1282 end subroutine af_consistent_fluxes
1283
1284 !> The neighbor nb has no children and id does, so set flux on the neighbor
1285 !> from our children. This ensures flux consistency at refinement boundary.
1286 subroutine flux_from_children(boxes, id, nb, f_ixs)
1287 type(box_t), intent(inout) :: boxes(:) !< List of all the boxes
1288 integer, intent(in) :: id !< Id of box for which we set fluxes
1289 integer, intent(in) :: nb !< Direction in which fluxes are set
1290 integer, intent(in) :: f_ixs(:) !< Indices of the fluxes
1291 integer :: nc, nch, c_id, i_ch, i, ic, d
1292 integer :: n_chnb, nb_id, i_nb
1293#if NDIM > 1
1294 integer :: ioff(NDIM)
1295#endif
1296#if NDIM == 2
1297 integer :: n
1298 real(dp) :: w1, w2
1299#endif
1300
1301
1302 nc = boxes(id)%n_cell
1303 nch = ishft(nc, -1) ! nc/2
1304 d = af_neighb_dim(nb)
1305 n_chnb = 2**(ndim-1)
1306 nb_id = boxes(id)%neighbors(nb)
1307
1308 if (af_neighb_low(nb)) then
1309 i = 1
1310 i_nb = nc+1
1311 else
1312 i = nc+1
1313 i_nb = 1
1314 end if
1315
1316 select case (d)
1317#if NDIM == 1
1318 case (1)
1319 do ic = 1, n_chnb
1320 ! Get index of child adjacent to neighbor
1321 i_ch = af_child_adj_nb(ic, nb)
1322 c_id = boxes(id)%children(i_ch)
1323 boxes(nb_id)%fc(i_nb, 1, f_ixs) = boxes(c_id)%fc(i, 1, f_ixs)
1324 end do
1325#elif NDIM == 2
1326 case (1)
1327 do ic = 1, n_chnb
1328 ! Get index of child adjacent to neighbor
1329 i_ch = af_child_adj_nb(ic, nb)
1330 c_id = boxes(id)%children(i_ch)
1331 ! Index offset of child w.r.t. parent
1332 ioff = nch*af_child_dix(:, i_ch)
1333 boxes(nb_id)%fc(i_nb, ioff(2)+1:ioff(2)+nch, 1, f_ixs) = 0.5_dp * ( &
1334 boxes(c_id)%fc(i, 1:nc:2, 1, f_ixs) + &
1335 boxes(c_id)%fc(i, 2:nc:2, 1, f_ixs))
1336 end do
1337 case (2)
1338 if (boxes(nb_id)%coord_t == af_cyl) then
1339 ! In cylindrical symmetry, we take the weighted average
1340 do ic = 1, n_chnb
1341 i_ch = af_child_adj_nb(ic, nb)
1342 c_id = boxes(id)%children(i_ch)
1343 ioff = nch*af_child_dix(:, i_ch)
1344
1345 do n = 1, nch
1346 call af_cyl_child_weights(boxes(nb_id), ioff(1)+n, w1, w2)
1347 boxes(nb_id)%fc(ioff(1)+n, i_nb, 2, f_ixs) = 0.5_dp * (&
1348 w1 * boxes(c_id)%fc(2*n-1, i, 2, f_ixs) + &
1349 w2 * boxes(c_id)%fc(2*n, i, 2, f_ixs))
1350 end do
1351 end do
1352 else
1353 ! Just take the average of the fine fluxes
1354 do ic = 1, n_chnb
1355 i_ch = af_child_adj_nb(ic, nb)
1356 c_id = boxes(id)%children(i_ch)
1357 ioff = nch*af_child_dix(:, i_ch)
1358 boxes(nb_id)%fc(ioff(1)+1:ioff(1)+nch, i_nb, 2, f_ixs) = 0.5_dp * ( &
1359 boxes(c_id)%fc(1:nc:2, i, 2, f_ixs) + &
1360 boxes(c_id)%fc(2:nc:2, i, 2, f_ixs))
1361 end do
1362 end if
1363#elif NDIM == 3
1364 case (1)
1365 do ic = 1, n_chnb
1366 i_ch = af_child_adj_nb(ic, nb)
1367 c_id = boxes(id)%children(i_ch)
1368 ioff = nch*af_child_dix(:, i_ch)
1369 boxes(nb_id)%fc(i_nb, ioff(2)+1:ioff(2)+nch, &
1370 ioff(3)+1:ioff(3)+nch, 1, f_ixs) = 0.25_dp * ( &
1371 boxes(c_id)%fc(i, 1:nc:2, 1:nc:2, 1, f_ixs) + &
1372 boxes(c_id)%fc(i, 2:nc:2, 1:nc:2, 1, f_ixs) + &
1373 boxes(c_id)%fc(i, 1:nc:2, 2:nc:2, 1, f_ixs) + &
1374 boxes(c_id)%fc(i, 2:nc:2, 2:nc:2, 1, f_ixs))
1375 end do
1376 case (2)
1377 do ic = 1, n_chnb
1378 i_ch = af_child_adj_nb(ic, nb)
1379 c_id = boxes(id)%children(i_ch)
1380 ioff = nch*af_child_dix(:, i_ch)
1381 boxes(nb_id)%fc(ioff(1)+1:ioff(1)+nch, i_nb, &
1382 ioff(3)+1:ioff(3)+nch, 2, f_ixs) = 0.25_dp * ( &
1383 boxes(c_id)%fc(1:nc:2, i, 1:nc:2, 2, f_ixs) + &
1384 boxes(c_id)%fc(2:nc:2, i, 1:nc:2, 2, f_ixs) + &
1385 boxes(c_id)%fc(1:nc:2, i, 2:nc:2, 2, f_ixs) + &
1386 boxes(c_id)%fc(2:nc:2, i, 2:nc:2, 2, f_ixs))
1387 end do
1388 case (3)
1389 do ic = 1, n_chnb
1390 i_ch = af_child_adj_nb(ic, nb)
1391 c_id = boxes(id)%children(i_ch)
1392 ioff = nch*af_child_dix(:, i_ch)
1393 boxes(nb_id)%fc(ioff(1)+1:ioff(1)+nch, &
1394 ioff(2)+1:ioff(2)+nch, i_nb, 3, f_ixs) = 0.25_dp * ( &
1395 boxes(c_id)%fc(1:nc:2, 1:nc:2, i, 3, f_ixs) + &
1396 boxes(c_id)%fc(2:nc:2, 1:nc:2, i, 3, f_ixs) + &
1397 boxes(c_id)%fc(1:nc:2, 2:nc:2, i, 3, f_ixs) + &
1398 boxes(c_id)%fc(2:nc:2, 2:nc:2, i, 3, f_ixs))
1399 end do
1400#endif
1401 end select
1402 end subroutine flux_from_children
1403
1404end module m_af_core
To fill ghost cells near physical boundaries in a custom way. If the number of ghost cells to fill is...
Definition: m_af_types.f90:477
To fill ghost cells near physical boundaries.
Definition: m_af_types.f90:464
To set cell-centered variables based on a user-defined function. This can be useful to avoid recomput...
Definition: m_af_types.f90:491
Subroutine for prolongation.
Definition: m_af_types.f90:498
To fill ghost cells near refinement boundaries.
Definition: m_af_types.f90:454
Subroutine for restriction.
Definition: m_af_types.f90:509
This module contains the core routines of Afivo, namely those that deal with initializing and changin...
Definition: m_af_core.f90:3
This module contains routines related to the filling of ghost cells. Note that corner ghost cells are...
Module containing slope limiters.
This module contains the routines related to prolongation: going from coarse to fine variables.
Definition: m_af_prolong.f90:3
This module contains routines for restriction: going from fine to coarse variables.
This module contains the basic types and constants that are used in the NDIM-dimensional version of A...
Definition: m_af_types.f90:3
This module contains all kinds of different 'helper' routines for Afivo. If the number of routines fo...
Definition: m_af_utils.f90:4
Type which stores all the boxes and levels, as well as some information about the number of boxes,...
Definition: m_af_types.f90:326
The basic building block of afivo: a box with cell-centered and face centered data,...
Definition: m_af_types.f90:286
Type which contains the indices of all boxes at a refinement level, as well as a list with all the "l...
Definition: m_af_types.f90:76