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 helmholtz_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 real(dp), parameter :: lambda = 1000.0_dp
19
20 type(af_t) :: tree
21 type(ref_info_t) :: ref_info
22 integer :: mg_iter
23 real(dp) :: residu(2), anal_err(2)
24 character(len=100) :: fname
25 type(mg_t) :: mg
26 type(gauss_t) :: gs
27 integer :: count_rate,t_start, t_end
28
29 print *, "Running helmholtz_cyl"
30 print *, "Number of threads", af_get_max_threads()
31
32
33 call gauss_init(gs, [1.0_dp, 1.0_dp], [0.01_dp, 0.04_dp], &
34 reshape([0.0_dp, 0.25_dp, 0.0_dp, 0.75_dp], [2,2]))
35
36 call af_add_cc_variable(tree, "phi", ix=i_phi)
37 call af_add_cc_variable(tree, "rhs", ix=i_rhs)
38 call af_add_cc_variable(tree, "err", ix=i_err)
39 call af_add_cc_variable(tree, "tmp", ix=i_tmp)
40
41
42 call af_init(tree, &
43 box_size, &
44 [1.0_dp, 1.0_dp], &
45 [box_size, box_size], &
46 coord=af_cyl)
47
48 call af_print_info(tree)
49
50 call system_clock(t_start, count_rate)
51 do
52
53 call af_loop_box(tree, set_init_cond)
54
55
56 call af_adjust_refinement(tree, ref_routine, ref_info)
57
58
59 if (ref_info%n_add == 0) exit
60 end do
61 call system_clock(t_end, count_rate)
62
63 write(*,"(A,Es10.3,A)") " Wall-clock time generating AMR grid: ", &
64 (t_end-t_start) / real(count_rate,dp), " seconds"
65
66 call af_print_info(tree)
67
68
69 mg%i_phi = i_phi
70 mg%i_rhs = i_rhs
71 mg%i_tmp = i_tmp
72 mg%sides_bc => sides_bc
73 mg%helmholtz_lambda = lambda
74
75
76
77
78
79 call mg_init(tree, mg)
80
81 print *, "Multigrid iteration | max residual | max error"
82 call system_clock(t_start, count_rate)
83
84 do mg_iter = 1, n_iterations
85
86
87
88 call mg_fas_fmg(tree, mg, .true., mg_iter>1)
89
90
91 call af_loop_box(tree, set_err)
92
93
94 call af_tree_min_cc(tree, i_tmp, residu(1))
95 call af_tree_max_cc(tree, i_tmp, residu(2))
96 call af_tree_min_cc(tree, i_err, anal_err(1))
97 call af_tree_max_cc(tree, i_err, anal_err(2))
98 write(*,"(I8,2Es14.5)") mg_iter, maxval(abs(residu)), &
99 maxval(abs(anal_err))
100
101 write(fname, "(A,I0)") "output/helmholtz_cyl_", mg_iter
102 call af_write_vtk(tree, trim(fname))
103 end do
104 call system_clock(t_end, count_rate)
105
106 write(*, "(A,I0,A,E10.3,A)") &
107 " Wall-clock time after ", n_iterations, &
108 " iterations: ", (t_end-t_start) / real(count_rate, dp), &
109 " seconds"
110
111
112
113 call af_destroy(tree)
114
115contains
116
117
118 subroutine ref_routine(box, cell_flags)
119 type(box_t), intent(in) :: box
120 integer, intent(out) :: cell_flags(box%n_cell, box%n_cell)
121 integer :: i, j, nc
122 real(dp) :: crv, dr2
123
124 nc = box%n_cell
125 dr2 = maxval(box%dr)**2
126
127
128 do j = 1, nc
129 do i = 1, nc
130 crv = dr2 * abs(box%cc(i, j, i_rhs))
131
132
133 if (crv > 5.0e-4_dp) then
134 cell_flags(i, j) = af_do_ref
135 else
136 cell_flags(i, j) = af_keep_ref
137 end if
138
139
140
141
142
143 end do
144 end do
145 end subroutine ref_routine
146
147
148 subroutine set_init_cond(box)
149 type(box_t), intent(inout) :: box
150 integer :: i, j, nc
151 real(dp) :: rz(2)
152
153 nc = box%n_cell
154
155 do j = 0, nc+1
156 do i = 0, nc+1
157 rz = af_r_cc(box, [i,j])
158 box%cc(i, j, i_rhs) = gauss_laplacian_cyl(gs, rz) - &
159 lambda * gauss_value(gs, rz)
160 end do
161 end do
162 end subroutine set_init_cond
163
164
165 subroutine set_err(box)
166 type(box_t), intent(inout) :: box
167 integer :: i, j, nc
168 real(dp) :: rz(2)
169
170 nc = box%n_cell
171 do j = 1, nc
172 do i = 1, nc
173 rz = af_r_cc(box, [i,j])
174 box%cc(i, j, i_err) = box%cc(i, j, i_phi) - gauss_value(gs, rz)
175 end do
176 end do
177 end subroutine set_err
178
179
180 subroutine sides_bc(box, nb, iv, coords, bc_val, bc_type)
181 type(box_t), intent(in) :: box
182 integer, intent(in) :: nb
183 integer, intent(in) :: iv
184 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
185 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
186 integer, intent(out) :: bc_type
187 integer :: n
188
189 if (nb == af_neighb_lowx) then
190
191 bc_type = af_bc_neumann
192 bc_val = 0.0_dp
193 else
194
195 bc_type = af_bc_dirichlet
196
197
198 do n = 1, box%n_cell**(ndim-1)
199 bc_val(n) = gauss_value(gs, coords(:, n))
200 end do
201 end if
202 end subroutine sides_bc
203
204end program helmholtz_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.