Afivo 0.3
All Classes Namespaces Functions Variables Pages
poisson_basic.f90

Example showing how to use multigrid and compare with an analytic solution. A standard 5-point Laplacian is used.

Example showing how to use multigrid and compare with an analytic solution. A standard 5-point Laplacian is used.

1#include "../src/cpp_macros.h"
2!> \example poisson_basic.f90
3!>
4!> Example showing how to use multigrid and compare with an analytic solution. A
5!> standard 5-point Laplacian is used.
6program poisson_basic
7 use m_af_all
9
10 implicit none
11
12 integer, parameter :: box_size = 16
13 integer, parameter :: n_iterations = 10
14 integer :: domain_size(NDIM)
15 real(dp) :: domain_len(NDIM)
16 integer :: i_phi
17 integer :: i_sol
18 integer :: i_rhs
19 integer :: i_err
20 integer :: i_tmp
21 integer :: i_gradx
22 integer :: i_egradx
23 type(af_t) :: tree
24 type(ref_info_t) :: refine_info
25 integer :: mg_iter
26 real(dp) :: residu(2), anal_err(2)
27 character(len=100) :: fname
28 type(mg_t) :: mg
29 type(gauss_t) :: gs
30 integer :: count_rate,t_start,t_end
31
32 print *, "Running poisson_basic_" // dimname
33 print *, "Number of threads", af_get_max_threads()
34
35 !> [Gauss_init]
36 ! The manufactured solution exists of Gaussians
37 ! Amplitudes: [1.0_dp, 1.0_dp]
38 ! Sigmas : [0.04_dp, 0.04_dp]
39 ! Locations : x, y, z = 0.25 or x, y, z = 0.75
40 call gauss_init(gs, [1.0_dp, 1.0_dp], [0.04_dp, 0.04_dp], &
41 reshape([dtimes(0.1_dp), &
42 dtimes(0.75_dp)], [ndim,2]))
43 !> [Gauss_init]
44
45 !> [af_init]
46 call af_add_cc_variable(tree, "phi", ix=i_phi)
47 call af_add_cc_variable(tree, "sol", ix=i_sol)
48 call af_add_cc_variable(tree, "rhs", ix=i_rhs)
49 call af_add_cc_variable(tree, "err", ix=i_err)
50 call af_add_cc_variable(tree, "tmp", ix=i_tmp)
51 call af_add_cc_variable(tree, "Dx", ix=i_gradx)
52 call af_add_cc_variable(tree, "eDx", ix=i_egradx)
53
54 domain_len(1) = 3.0_dp
55 domain_len(2:) = 1.0_dp
56 domain_size(1) = 3 * box_size
57 domain_size(2:) = box_size
58
59 ! Initialize tree
60 call af_init(tree, & ! Tree to initialize
61 box_size, & ! A box contains box_size**DIM cells
62 domain_len, &
63 domain_size)
64 !> [af_init]
65
66 call system_clock(t_start, count_rate)
67 !> [set_refinement]
68 do
69 ! For each box, set the initial conditions
70 call af_loop_box(tree, set_initial_condition)
71
72 ! This updates the refinement of the tree, by at most one level per call.
73 ! The second argument is a subroutine that is called for each box that can
74 ! be refined or derefined, and it should set refinement flags. Information
75 ! about the changes in refinement are returned in the third argument.
76 call af_adjust_refinement(tree, refine_routine, refine_info)
77
78 ! If no new boxes have been added, exit the loop
79 if (refine_info%n_add == 0) exit
80 end do
81 !> [set_refinement]
82 call system_clock(t_end, count_rate)
83
84 write(*,"(A,Es10.3,A)") " Wall-clock time generating AMR grid: ", &
85 (t_end-t_start) / real(count_rate,dp), " seconds"
86
87 call af_print_info(tree)
88
89 ! Initialize the multigrid options. This performs some basics checks and sets
90 ! default values where necessary.
91 mg%i_phi = i_phi ! Solution variable
92 mg%i_rhs = i_rhs ! Right-hand side variable
93 mg%i_tmp = i_tmp ! Variable for temporary space
94 mg%sides_bc => sides_bc ! Method for boundary conditions
95
96 ! This routine does not initialize the multigrid fields boxes%i_phi,
97 ! boxes%i_rhs and boxes%i_tmp. These fileds will be initialized at the
98 ! first call of mg_fas_fmg
99 call mg_init(tree, mg)
100
101 print *, "Multigrid iteration | max residual | max error"
102 call system_clock(t_start, count_rate)
103
104 do mg_iter = 1, n_iterations
105 ! Perform a FAS-FMG cycle (full approximation scheme, full multigrid). The
106 ! third argument controls whether the residual is stored in i_tmp. The
107 ! fourth argument controls whether to improve the current solution.
108 call mg_fas_fmg(tree, mg, set_residual=.true., have_guess=(mg_iter>1))
109
110 ! Compute the error compared to the analytic solution
111 call af_loop_box(tree, set_error)
112
113 ! Determine the minimum and maximum residual and error
114 call af_tree_min_cc(tree, i_tmp, residu(1))
115 call af_tree_max_cc(tree, i_tmp, residu(2))
116 call af_tree_min_cc(tree, i_err, anal_err(1))
117 call af_tree_max_cc(tree, i_err, anal_err(2))
118 write(*,"(I8,13x,2(Es14.5))") mg_iter, maxval(abs(residu)), &
119 maxval(abs(anal_err))
120
121 !> [write_output]
122 ! This writes a Silo output file containing the cell-centered values of the
123 ! leaves of the tree (the boxes not covered by refinement).
124 write(fname, "(A,I0)") "output/poisson_basic_" // dimname // "_", mg_iter
125 call af_write_silo(tree, trim(fname))
126 !> [write_output]
127 end do
128 call system_clock(t_end, count_rate)
129
130 write(*, "(A,I0,A,E10.3,A)") &
131 " Wall-clock time after ", n_iterations, &
132 " iterations: ", (t_end-t_start) / real(count_rate, dp), &
133 " seconds"
134
135 ! This call is not really necessary here, but cleaning up the data in a tree
136 ! is important if your program continues with other tasks.
137 call af_destroy(tree)
138
139contains
140
141 !> [refine_routine]
142 ! Return the refinement flags for box
143 subroutine refine_routine(box, cell_flags)
144 type(box_t), intent(in) :: box
145 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
146 integer :: IJK, nc
147 real(dp) :: rr(NDIM), dr2, drhs
148
149 nc = box%n_cell
150 dr2 = maxval(box%dr)**2
151
152 do kji_do(1,nc)
153 rr = af_r_cc(box, [ijk])
154
155 ! This is an estimate of the truncation error in the right-hand side,
156 ! which is related to the fourth derivative of the solution.
157 drhs = dr2 * box%cc(ijk, i_rhs)
158
159 if (abs(drhs) > 1e-3_dp .and. box%lvl < 7 - ndim) then
160 cell_flags(ijk) = af_do_ref
161 else
162 cell_flags(ijk) = af_keep_ref
163 end if
164 end do; close_do
165 end subroutine refine_routine
166 !> [refine_routine]
167
168 !> [set_initial_condition]
169 ! This routine sets the initial conditions for each box
170 subroutine set_initial_condition(box)
171 type(box_t), intent(inout) :: box
172 integer :: IJK, nc
173 real(dp) :: rr(NDIM), grad(NDIM)
174
175 nc = box%n_cell
176
177 do kji_do(1,nc)
178 ! Get the coordinate of the cell center at IJK
179 rr = af_r_cc(box, [ijk])
180
181 ! And set the rhs values
182 box%cc(ijk, i_rhs) = gauss_laplacian(gs, rr)
183 call gauss_gradient(gs, rr, grad)
184 box%cc(ijk, i_gradx) = grad(1)
185 box%cc(ijk, i_sol) = gauss_value(gs, rr)
186 end do; close_do
187 end subroutine set_initial_condition
188 !> [set_initial_condition]
189
190 !> [set_error]
191 ! Set the error compared to the analytic solution
192 subroutine set_error(box)
193 type(box_t), intent(inout) :: box
194 integer :: IJK, nc
195 real(dp) :: rr(NDIM), gradx
196
197 nc = box%n_cell
198
199 do kji_do(1,nc)
200 rr = af_r_cc(box, [ijk])
201 box%cc(ijk, i_err) = box%cc(ijk, i_phi) - gauss_value(gs, rr)
202#if NDIM == 1
203 gradx = 0.5_dp * (box%cc(i+1, i_phi) - &
204 box%cc(i-1, i_phi)) / box%dr(1)
205#elif NDIM == 2
206 gradx = 0.5_dp * (box%cc(i+1, j, i_phi) - &
207 box%cc(i-1, j, i_phi)) / box%dr(1)
208#elif NDIM == 3
209 gradx = 0.5_dp * (box%cc(i+1, j, k, i_phi) - &
210 box%cc(i-1, j, k, i_phi)) / box%dr(1)
211#endif
212 box%cc(ijk, i_egradx) = gradx - box%cc(ijk, i_gradx)
213
214 end do; close_do
215 end subroutine set_error
216 !> [set_error]
217
218 ! This routine sets boundary conditions for a box
219 subroutine sides_bc(box, nb, iv, coords, bc_val, bc_type)
220 type(box_t), intent(in) :: box
221 integer, intent(in) :: nb
222 integer, intent(in) :: iv
223 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
224 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
225 integer, intent(out) :: bc_type
226 integer :: n
227
228 ! We use dirichlet boundary conditions
229 bc_type = af_bc_dirichlet
230
231 ! Below the solution is specified in the approriate ghost cells
232 do n = 1, box%n_cell**(ndim-1)
233 bc_val(n) = gauss_value(gs, coords(:, n))
234 end do
235 end subroutine sides_bc
236
237end program poisson_basic
Module which contains all Afivo modules, so that a user does not have to include them separately.
Definition: m_af_all.f90:3
This module can be used to construct solutions consisting of one or more Gaussians.
Definition: m_gaussians.f90:3