Example showing how to use multigrid with a numerical stencil.
1#include "../src/cpp_macros.h"
2
3
4
5program poisson_stencil
8
9 implicit none
10
11 integer, parameter :: box_size = 16
12 integer, parameter :: n_iterations = 10
13 integer :: domain_size(NDIM)
14 real(dp) :: domain_len(NDIM)
15 integer :: i_phi
16 integer :: i_sol
17 integer :: i_rhs
18 integer :: i_err
19 integer :: i_tmp
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_stencil_" // 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.1_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, "sol", ix=i_sol)
38 call af_add_cc_variable(tree, "rhs", ix=i_rhs)
39 call af_add_cc_variable(tree, "err", ix=i_err)
40 call af_add_cc_variable(tree, "tmp", ix=i_tmp)
41
42 domain_len(1) = 3.0_dp
43 domain_len(2:) = 1.0_dp
44 domain_size(1) = 3 * box_size
45 domain_size(2:) = box_size
46
47 call af_init(tree, &
48 box_size, &
49 domain_len, &
50 domain_size)
51
52 call system_clock(t_start, count_rate)
53 do
54 call af_loop_box(tree, set_initial_condition)
55 call af_adjust_refinement(tree, refine_routine, refine_info)
56 if (refine_info%n_add == 0) exit
57 end do
58 call system_clock(t_end, count_rate)
59
60 write(*,"(A,Es10.3,A)") " Wall-clock time generating AMR grid: ", &
61 (t_end-t_start) / real(count_rate,dp), " seconds"
62
63 call af_print_info(tree)
64
65 mg%i_phi = i_phi
66 mg%i_rhs = i_rhs
67 mg%i_tmp = i_tmp
68 mg%sides_bc => sides_bc
69
70 mg%operator_key = mg_normal_box
71
72 call mg_init(tree, mg)
73
74 print *, "Multigrid iteration | max residual | max error"
75 call system_clock(t_start, count_rate)
76
77 do mg_iter = 1, n_iterations
78 call mg_fas_fmg(tree, mg, set_residual=.true., have_guess=(mg_iter>1))
79
80 call af_loop_box(tree, set_error)
81
82 call af_tree_min_cc(tree, i_tmp, residu(1))
83 call af_tree_max_cc(tree, i_tmp, residu(2))
84 call af_tree_min_cc(tree, i_err, anal_err(1))
85 call af_tree_max_cc(tree, i_err, anal_err(2))
86 write(*,"(I8,13x,2(Es14.5))") mg_iter, maxval(abs(residu)), &
87 maxval(abs(anal_err))
88
89 write(fname, "(A,I0)") "output/poisson_stencil_" // dimname // "_", mg_iter
90 call af_write_silo(tree, trim(fname))
91 end do
92 call system_clock(t_end, count_rate)
93
94 write(*, "(A,I0,A,E10.3,A)") &
95 " Wall-clock time after ", n_iterations, &
96 " iterations: ", (t_end-t_start) / real(count_rate, dp), &
97 " seconds"
98
99 call af_destroy(tree)
100
101contains
102
103 subroutine refine_routine(box, cell_flags)
104 type(box_t), intent(in) :: box
105 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
106 integer :: IJK, nc
107 real(dp) :: rr(NDIM), dr2, drhs
108
109 nc = box%n_cell
110 dr2 = maxval(box%dr)**2
111
112 do kji_do(1,nc)
113 rr = af_r_cc(box, [ijk])
114
115 drhs = dr2 * box%cc(ijk, i_rhs)
116
117 if (abs(drhs) > 1e-3_dp .and. box%lvl < 7 - ndim) then
118 cell_flags(ijk) = af_do_ref
119 else
120 cell_flags(ijk) = af_keep_ref
121 end if
122 end do; close_do
123 end subroutine refine_routine
124
125 subroutine set_initial_condition(box)
126 type(box_t), intent(inout) :: box
127 integer :: IJK, nc
128 real(dp) :: rr(NDIM)
129
130 nc = box%n_cell
131 do kji_do(1,nc)
132 rr = af_r_cc(box, [ijk])
133 box%cc(ijk, i_rhs) = gauss_laplacian(gs, rr)
134 box%cc(ijk, i_sol) = gauss_value(gs, rr)
135 end do; close_do
136 end subroutine set_initial_condition
137
138 subroutine set_error(box)
139 type(box_t), intent(inout) :: box
140 integer :: IJK, nc
141 real(dp) :: rr(NDIM)
142
143 nc = box%n_cell
144
145 do kji_do(1,nc)
146 rr = af_r_cc(box, [ijk])
147 box%cc(ijk, i_err) = box%cc(ijk, i_phi) - gauss_value(gs, rr)
148 end do; close_do
149 end subroutine set_error
150
151 subroutine sides_bc(box, nb, iv, coords, bc_val, bc_type)
152 type(box_t), intent(in) :: box
153 integer, intent(in) :: nb
154 integer, intent(in) :: iv
155 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
156 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
157 integer, intent(out) :: bc_type
158 integer :: n
159
160 bc_type = af_bc_dirichlet
161
162 do n = 1, box%n_cell**(ndim-1)
163 bc_val(n) = gauss_value(gs, coords(:, n))
164 end do
165 end subroutine sides_bc
166
167end program poisson_stencil
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.