Example showing how to use multigrid in combination with Neumann boundary conditions. There is also a test case for a cylindrical geometry.
Example showing how to use multigrid in combination with Neumann boundary conditions. There is also a test case for a cylindrical geometry.
1#include "../src/cpp_macros.h"
2
3
4
5
6program poisson_neumann
8
9 implicit none
10
11 integer, parameter :: box_size = 8
12 integer, parameter :: n_iterations = 10
13 integer :: i_err
14
15 type(af_t) :: tree
16 type(ref_info_t) :: refine_info
17 integer :: mg_iter
18 character(len=100) :: fname, arg
19 type(mg_t) :: mg
20 integer :: count_rate, t_start, t_end, coord
21 logical :: cylindrical = .false.
22
23 print *, "Running poisson_neumann_" // dimname // ""
24 print *, "Number of threads", af_get_max_threads()
25
26 if (command_argument_count() > 0) then
27 call get_command_argument(1, arg)
28 read (arg, *) cylindrical
29 end if
30
31 if (cylindrical) then
32 coord = af_cyl
33 else
34 coord = af_xyz
35 end if
36
37 call af_add_cc_variable(tree, "phi", ix=mg%i_phi)
38 call af_add_cc_variable(tree, "rhs", ix=mg%i_rhs)
39 call af_add_cc_variable(tree, "tmp", ix=mg%i_tmp)
40 call af_add_cc_variable(tree, "err", ix=i_err)
41
42
43 call af_init(tree, &
44 box_size, &
45 [dtimes(1.0_dp)], &
46 [dtimes(box_size)], &
47 coord=coord)
48
49 call system_clock(t_start, count_rate)
50 do
51 call af_loop_box(tree, set_initial_condition)
52 call af_adjust_refinement(tree, refine_routine, refine_info, 0)
53 if (refine_info%n_add == 0) exit
54 end do
55 call system_clock(t_end, count_rate)
56
57 write(*,"(A,Es10.3,A)") " Wall-clock time generating AMR grid: ", &
58 (t_end-t_start) / real(count_rate,dp), " seconds"
59
60 call af_print_info(tree)
61
62
63 if (cylindrical) then
64 mg%sides_bc => sides_bc_cyl
65 else
66 mg%sides_bc => sides_bc
67 end if
68
69 call mg_init(tree, mg)
70
71 print *, "Multigrid iteration | max residual | max error"
72 call system_clock(t_start, count_rate)
73
74 do mg_iter = 1, n_iterations
75 call mg_fas_fmg(tree, mg, set_residual=.true., have_guess=(mg_iter>1))
76
77 call af_loop_box(tree, set_error)
78
79 write(fname, "(A,I0)") "output/poisson_neumann_" // dimname // "_", mg_iter
80 call af_write_silo(tree, trim(fname))
81 end do
82 call system_clock(t_end, count_rate)
83
84 write(*, "(A,I0,A,E10.3,A)") &
85 " Wall-clock time after ", n_iterations, &
86 " iterations: ", (t_end-t_start) / real(count_rate, dp), &
87 " seconds"
88
89contains
90
91
92 subroutine refine_routine(box, cell_flags)
93 type(box_t), intent(in) :: box
94 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
95
96
97 if (box%lvl <= 4 .and. all(box%r_min < 0.25_dp)) then
98 cell_flags(dtimes(:)) = af_do_ref
99 else
100 cell_flags(dtimes(:)) = af_keep_ref
101 end if
102 end subroutine refine_routine
103
104
105 subroutine set_initial_condition(box)
106 type(box_t), intent(inout) :: box
107
108 if (cylindrical) then
109
110
111 box%cc(dtimes(:), mg%i_rhs) = 4.0_dp
112 else
113 box%cc(dtimes(:), mg%i_rhs) = 0.0_dp
114 end if
115 end subroutine set_initial_condition
116
117
118 subroutine set_error(box)
119 type(box_t), intent(inout) :: box
120 integer :: IJK, nc
121 real(dp) :: rr(NDIM)
122
123 nc = box%n_cell
124 do kji_do(1,nc)
125 rr = af_r_cc(box, [ijk])
126 box%cc(ijk, i_err) = box%cc(ijk, mg%i_phi) - solution(rr(1))
127 end do; close_do
128 end subroutine set_error
129
130 real(dp) elemental function solution(x)
131 real(dp), intent(in) :: x
132
133 if (cylindrical) then
134 solution = x**2
135 else
136 solution = x
137 end if
138 end function solution
139
140
141 subroutine sides_bc(box, nb, iv, coords, bc_val, bc_type)
142 type(box_t), intent(in) :: box
143 integer, intent(in) :: nb
144 integer, intent(in) :: iv
145 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
146 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
147 integer, intent(out) :: bc_type
148 integer :: nc
149
150 nc = box%n_cell
151
152
153 select case (nb)
154 case (af_neighb_lowx)
155 bc_type = af_bc_dirichlet
156 bc_val = 0.0_dp
157 case (af_neighb_highx)
158 bc_type = af_bc_neumann
159 bc_val = 1.0_dp
160 case default
161 bc_type = af_bc_neumann
162 bc_val = 0.0_dp
163 end select
164 end subroutine sides_bc
165
166
167 subroutine sides_bc_cyl(box, nb, iv, coords, bc_val, bc_type)
168 type(box_t), intent(in) :: box
169 integer, intent(in) :: nb
170 integer, intent(in) :: iv
171 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
172 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
173 integer, intent(out) :: bc_type
174 integer :: nc
175
176 nc = box%n_cell
177
178
179 select case (nb)
180 case (af_neighb_lowx)
181 bc_type = af_bc_neumann
182 bc_val = 0.0_dp
183 case (af_neighb_highx)
184 bc_type = af_bc_neumann
185 bc_val = 2.0_dp
186 case default
187
188 bc_type = af_bc_dirichlet
189 bc_val = solution(coords(1, :))
190 end select
191 end subroutine sides_bc_cyl
192
193end program
Module which contains all Afivo modules, so that a user does not have to include them separately.