Example showing how to use multigrid and compare with an analytic solution, using the method of manufactured solutions. A standard 5-point Laplacian is used in cylindrical coordinates.
Example showing how to use multigrid and compare with an analytic solution, using the method of manufactured solutions. A standard 5-point Laplacian is used in cylindrical coordinates.
1
2
3
4
5
6program poisson_cyl
9
10 implicit none
11
12 integer, parameter :: box_size = 8
13 integer, parameter :: n_iterations = 10
14 integer :: i_phi
15 integer :: i_rhs
16 integer :: i_err
17 integer :: i_tmp
18
19 type(af_t) :: tree
20 type(ref_info_t) :: ref_info
21 integer :: mg_iter
22 real(dp) :: residu(2), anal_err(2)
23 character(len=100) :: fname
24 type(mg_t) :: mg
25 type(gauss_t) :: gs
26 integer :: count_rate,t_start, t_end
27
28 print *, "Running poisson_cyl"
29 print *, "Number of threads", af_get_max_threads()
30
31
32 call gauss_init(gs, [1.0_dp], [0.2_dp], &
33 reshape([0.0_dp, 0.5_dp], [2,1]))
34
35 call af_add_cc_variable(tree, "phi", ix=i_phi)
36 call af_add_cc_variable(tree, "rhs", ix=i_rhs)
37 call af_add_cc_variable(tree, "err", ix=i_err)
38 call af_add_cc_variable(tree, "tmp", ix=i_tmp)
39
40
41 call af_init(tree, &
42 box_size, &
43 [1.0_dp, 1.0_dp], &
44 [box_size, box_size], &
45 coord=af_cyl)
46
47 call af_print_info(tree)
48
49 call system_clock(t_start, count_rate)
50 do
51
52 call af_loop_box(tree, set_init_cond)
53
54
55 call af_adjust_refinement(tree, ref_routine, ref_info)
56
57
58 if (ref_info%n_add == 0) exit
59 end do
60 call system_clock(t_end, count_rate)
61
62 write(*,"(A,Es10.3,A)") " Wall-clock time generating AMR grid: ", &
63 (t_end-t_start) / real(count_rate,dp), " seconds"
64
65 call af_print_info(tree)
66
67
68 mg%i_phi = i_phi
69 mg%i_rhs = i_rhs
70 mg%i_tmp = i_tmp
71 mg%sides_bc => sides_bc
72
73
74
75
76
77 call mg_init(tree, mg)
78
79 print *, "Multigrid iteration | max residual | max error"
80 call system_clock(t_start, count_rate)
81
82 do mg_iter = 1, n_iterations
83
84
85
86 call mg_fas_fmg(tree, mg, .true., mg_iter>1)
87
88
89 call af_loop_box(tree, set_err)
90
91
92 call af_tree_min_cc(tree, i_tmp, residu(1))
93 call af_tree_max_cc(tree, i_tmp, residu(2))
94 call af_tree_min_cc(tree, i_err, anal_err(1))
95 call af_tree_max_cc(tree, i_err, anal_err(2))
96 write(*,"(I8,2Es14.5)") mg_iter, maxval(abs(residu)), &
97 maxval(abs(anal_err))
98
99 write(fname, "(A,I0)") "output/poisson_cyl_", mg_iter
100 call af_write_silo(tree, trim(fname))
101 end do
102 call system_clock(t_end, count_rate)
103
104 write(*, "(A,I0,A,E10.3,A)") &
105 " Wall-clock time after ", n_iterations, &
106 " iterations: ", (t_end-t_start) / real(count_rate, dp), &
107 " seconds"
108
109
110
111 call af_destroy(tree)
112 call mg_destroy(mg)
113
114contains
115
116
117 subroutine ref_routine(box, cell_flags)
118 type(box_t), intent(in) :: box
119 integer, intent(out) :: cell_flags(box%n_cell, box%n_cell)
120 integer :: i, j, nc
121 real(dp) :: crv, dr2
122
123 nc = box%n_cell
124 dr2 = maxval(box%dr)**2
125
126
127 do j = 1, nc
128 do i = 1, nc
129 crv = dr2 * abs(box%cc(i, j, i_rhs))
130
131
132 if (crv > 5.0e-4_dp) then
133 cell_flags(i, j) = af_do_ref
134 else
135 cell_flags(i, j) = af_keep_ref
136 end if
137 end do
138 end do
139 end subroutine ref_routine
140
141
142 subroutine set_init_cond(box)
143 type(box_t), intent(inout) :: box
144 integer :: i, j, nc
145 real(dp) :: rz(2)
146
147 nc = box%n_cell
148
149 do j = 0, nc+1
150 do i = 0, nc+1
151 rz = af_r_cc(box, [i,j])
152 box%cc(i, j, i_rhs) = gauss_laplacian_cyl(gs, rz)
153 end do
154 end do
155 end subroutine set_init_cond
156
157
158 subroutine set_err(box)
159 type(box_t), intent(inout) :: box
160 integer :: i, j, nc
161 real(dp) :: rz(2)
162
163 nc = box%n_cell
164 do j = 1, nc
165 do i = 1, nc
166 rz = af_r_cc(box, [i,j])
167 box%cc(i, j, i_err) = box%cc(i, j, i_phi) - gauss_value(gs, rz)
168 end do
169 end do
170 end subroutine set_err
171
172
173 subroutine sides_bc(box, nb, iv, coords, bc_val, bc_type)
174 type(box_t), intent(in) :: box
175 integer, intent(in) :: nb
176 integer, intent(in) :: iv
177 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
178 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
179 integer, intent(out) :: bc_type
180 integer :: n
181
182 if (nb == af_neighb_lowx) then
183
184 bc_type = af_bc_neumann
185 bc_val = 0.0_dp
186 else
187
188 bc_type = af_bc_dirichlet
189
190
191 do n = 1, box%n_cell**(ndim-1)
192 bc_val(n) = gauss_value(gs, coords(:, n))
193 end do
194 end if
195 end subroutine sides_bc
196
197end program poisson_cyl
Module which contains all Afivo modules, so that a user does not have to include them separately.
This module can be used to construct solutions consisting of one or more Gaussians.