This example shows how to create an AMR tree, perform random refinement, write output files, and how to fill ghost cells.
This example shows how to create an AMR tree, perform random refinement, write output files, and how to fill ghost cells.
1#include "../src/cpp_macros.h"
2
3
4
5
6program random_refinement
8
9 implicit none
10
11 type(af_t) :: tree
12 integer :: grid_size(NDIM)
13 real(dp) :: domain_size(NDIM)
14 logical :: periodic(NDIM) = .true.
15 integer :: iter, boxes_used
16 integer, parameter :: coord_type = af_xyz
17
18 integer, parameter :: box_size = 8
19 integer, parameter :: i_phi = 1
20 type(ref_info_t) :: ref_info
21 real(dp) :: sum_phi_t0, sum_phi
22 character(len=100) :: fname
23 integer :: count_rate,t_start,t_end
24
25 write(*,'(A,I0,A)') 'program random_refinement_', ndim, "d"
26
27 print *, "Number of threads", af_get_max_threads()
28
29 call af_add_cc_variable(tree, "phi")
30
31 call af_set_cc_methods(tree, 1, af_bc_dirichlet_zero, &
32 prolong=af_prolong_linear)
33
34
35 grid_size(1) = 2 * box_size
36 grid_size(2:) = box_size
37 domain_size(1) = 2 * acos(-1.0_dp)
38 domain_size(2:) = acos(-1.0_dp)
39
40 if (command_argument_count() == 1) then
41
42 call get_command_argument(1, fname)
43 call af_read_tree(tree, fname, read_other_data=read_string)
44 else
45 call af_init(tree, &
46 box_size, &
47 domain_size, &
48 grid_size, &
49 periodic=periodic, &
50 coord=coord_type)
51
52 call af_print_info(tree)
53
54
55
56
57 call af_loop_box(tree, set_init_cond)
58
59
60
61
62 call af_gc_tree(tree, [i_phi])
63 end if
64
65 call af_tree_sum_cc(tree, i_phi, sum_phi_t0)
66
67 call system_clock(t_start, count_rate)
68 boxes_used = 1
69 do iter = 1, 50
70
71
72
73 write(fname, "(A,I0)") "output/random_refinement_" // dimname // "_", iter
74 call af_write_silo(tree, trim(fname), n_cycle=iter)
75 call af_write_tree(tree, fname, write_other_data=write_string)
76
77
78
79
80
81
82
83
84
85
86 call af_adjust_refinement(tree, ref_routine, ref_info, ref_buffer=0)
87
88 call af_tree_sum_cc(tree, i_phi, sum_phi)
89 write(*, "(A,E10.2)") " conservation error: ", sum_phi - sum_phi_t0
90 boxes_used = boxes_used + ref_info%n_add - ref_info%n_rm
91 write(*,'(4(3x,A,1x,i6))') "# new boxes", ref_info%n_add, &
92 "# removed boxes", ref_info%n_rm, &
93 "# boxes used ", boxes_used, &
94 " highest level ", tree%highest_lvl
95 print *, "boxes / in use / not used: ", tree%highest_id, &
96 count(tree%boxes(1:tree%highest_id)%in_use), &
97 tree%n_removed_ids
98 end do
99 call system_clock(t_end, count_rate)
100
101 write(*, '(A,i3,1x,A,f8.2,1x,A,/)') &
102 ' Wall-clock time after ',iter, &
103 ' iterations: ', (t_end-t_start) / real(count_rate, dp), &
104 ' seconds'
105
106 call af_print_info(tree)
107
108
109
110 call af_destroy(tree)
111
112contains
113
114
115 subroutine ref_routine(box, cell_flags)
116 type(box_t), intent(in) :: box
117 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
118 real(dp) :: rr
119
120
121 call random_number(rr)
122
123 if (rr < 0.5_dp**ndim .and. box%lvl < 5) then
124 cell_flags = af_do_ref
125 else
126 cell_flags = af_rm_ref
127
128 end if
129 end subroutine ref_routine
130
131
132 subroutine set_init_cond(box)
133 type(box_t), intent(inout) :: box
134 integer :: IJK, nc
135 real(dp) :: rr(NDIM)
136
137 nc = box%n_cell
138 do kji_do(1,nc)
139
140 rr = af_r_cc(box, [ijk])
141
142
143#if NDIM > 1
144 box%cc(ijk, i_phi) = sin(0.5_dp * rr(1))**2 * cos(rr(2))**2
145#else
146 box%cc(ijk, i_phi) = sin(0.5_dp * rr(1))**2
147#endif
148 end do; close_do
149 end subroutine set_init_cond
150
151 subroutine write_string(my_unit)
152 integer, intent(in) :: my_unit
153 character(len=10) :: str
154
155 str = "Hello"
156 write(my_unit) str
157 end subroutine write_string
158
159 subroutine read_string(my_unit)
160 integer, intent(in) :: my_unit
161 character(len=10) :: str
162 read(my_unit) str
163
164 if (str /= "Hello") error stop "error reading other data"
165 end subroutine read_string
166
167end program random_refinement
Module which contains all Afivo modules, so that a user does not have to include them separately.