Example showing how to use multigrid for Helmholtz equation and compare with an analytic solution. A standard 5 or 7 point stencil is used.
Example showing how to use multigrid for Helmholtz equation and compare with an analytic solution. A standard 5 or 7 point stencil is used.
1#include "../src/cpp_macros.h"
2
3
4
5
6program poisson_helmholtz
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 = 1.0e3_dp
19
20 type(af_t) :: tree
21 type(ref_info_t) :: refine_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 poisson_helmholtz_" // dimname // ""
30 print *, "Number of threads", af_get_max_threads()
31
32 call gauss_init(gs, [1.0_dp, 1.0_dp], [0.04_dp, 0.04_dp], &
33 reshape([dtimes(0.25_dp), &
34 dtimes(0.75_dp)], [ndim,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, "tmp", ix=i_tmp)
39 call af_add_cc_variable(tree, "err", ix=i_err)
40
41
42 call af_init(tree, &
43 box_size, &
44 [dtimes(1.0_dp)], &
45 [dtimes(box_size)])
46
47 do
48 call af_loop_box(tree, set_initial_condition)
49 call af_adjust_refinement(tree, refine_routine, refine_info)
50 if (refine_info%n_add == 0) exit
51 end do
52
53 call af_print_info(tree)
54
55 mg%i_phi = i_phi
56 mg%i_rhs = i_rhs
57 mg%i_tmp = i_tmp
58 mg%sides_bc => sides_bc
59 mg%helmholtz_lambda = lambda
60
61 call mg_init(tree, mg)
62
63 print *, "Multigrid iteration | max residual | max error"
64 call system_clock(t_start, count_rate)
65
66 do mg_iter = 1, n_iterations
67 call mg_fas_fmg(tree, mg, set_residual=.true., have_guess=(mg_iter>1))
68 call af_loop_box(tree, set_error)
69 call af_tree_min_cc(tree, i_tmp, residu(1))
70 call af_tree_max_cc(tree, i_tmp, residu(2))
71 call af_tree_min_cc(tree, i_err, anal_err(1))
72 call af_tree_max_cc(tree, i_err, anal_err(2))
73 write(*,"(I8,13x,2(Es14.5))") mg_iter, maxval(abs(residu)), &
74 maxval(abs(anal_err))
75
76 write(fname, "(A,I0)") "output/poisson_helmholtz_" // dimname // "_", mg_iter
77 call af_write_vtk(tree, trim(fname))
78 end do
79 call system_clock(t_end, count_rate)
80
81 write(*, "(A,I0,A,E10.3,A)") &
82 " Wall-clock time after ", n_iterations, &
83 " iterations: ", (t_end-t_start) / real(count_rate, dp), &
84 " seconds"
85
86contains
87
88
89 subroutine refine_routine(box, cell_flags)
90 type(box_t), intent(in) :: box
91 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
92 integer :: IJK, nc
93 real(dp) :: rr(NDIM), dr2, drhs
94
95 nc = box%n_cell
96 dr2 = maxval(box%dr)**2
97
98 do kji_do(1,nc)
99 rr = af_r_cc(box, [ijk])
100
101
102
103 drhs = dr2 * box%cc(ijk, i_rhs)
104
105 if (abs(drhs) > 1e-3_dp .and. box%lvl < 5) then
106 cell_flags(ijk) = af_do_ref
107 else
108 cell_flags(ijk) = af_keep_ref
109 end if
110 end do; close_do
111 end subroutine refine_routine
112
113
114 subroutine set_initial_condition(box)
115 type(box_t), intent(inout) :: box
116 integer :: IJK, nc
117 real(dp) :: rr(NDIM)
118
119 nc = box%n_cell
120
121 do kji_do(1,nc)
122
123 rr = af_r_cc(box, [ijk])
124
125
126 box%cc(ijk, i_rhs) = gauss_laplacian(gs, rr) - &
127 lambda * gauss_value(gs, rr)
128 end do; close_do
129 end subroutine set_initial_condition
130
131
132 subroutine set_error(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 rr = af_r_cc(box, [ijk])
140 box%cc(ijk, i_err) = box%cc(ijk, i_phi) - gauss_value(gs, rr)
141 end do; close_do
142 end subroutine set_error
143
144
145 subroutine sides_bc(box, nb, iv, coords, bc_val, bc_type)
146 type(box_t), intent(in) :: box
147 integer, intent(in) :: nb
148 integer, intent(in) :: iv
149 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
150 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
151 integer, intent(out) :: bc_type
152 integer :: n
153
154
155 bc_type = af_bc_dirichlet
156
157
158 do n = 1, box%n_cell**(ndim-1)
159 bc_val(n) = gauss_value(gs, coords(:, n))
160 end do
161 end subroutine sides_bc
162
163end program poisson_helmholtz
164
165
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.