Afivo 0.3
All Classes Namespaces Functions Variables Pages
poisson_cyl_analytic.f90

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

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

1!> \example poisson_cyl_analytic.f90
2!>
3!> Example showing how to use multigrid and compare with an analytic solution. A
4!> standard 5-point Laplacian is used in cylindrical coordinates.
5program poisson_cyl_analytic
6 use m_af_all
7
8 implicit none
9
10 integer, parameter :: box_size = 8
11 integer, parameter :: n_iterations = 10
12 integer :: i_phi
13 integer :: i_rhs
14 integer :: i_err
15 integer :: i_sol
16 integer :: i_tmp
17 integer :: i_gradx
18 integer :: i_egradx
19
20 real(dp), parameter :: domain_len = 1.25e-2_dp
21 real(dp), parameter :: pi = acos(-1.0_dp)
22 real(dp), parameter :: sigma = 4e-4_dp * sqrt(0.5_dp)
23 real(dp), parameter :: rz_source(2) = [0.0_dp, 0.5_dp] * domain_len
24 real(dp), parameter :: epsilon_source = 8.85e-12_dp
25 real(dp), parameter :: Q_source = 3e18_dp * 1.6022e-19_dp * &
26 sigma**3 * sqrt(2 * pi)**3
27
28 type(af_t) :: tree
29 type(ref_info_t) :: ref_info
30 integer :: mg_iter
31 real(dp) :: residu(2), anal_err(2), max_val
32 character(len=100) :: fname
33 type(mg_t) :: mg
34 integer :: count_rate,t_start, t_end
35
36 print *, "Running poisson_cyl_analytic"
37 print *, "Number of threads", af_get_max_threads()
38
39 call af_add_cc_variable(tree, "phi", ix=i_phi)
40 call af_add_cc_variable(tree, "rhs", ix=i_rhs)
41 call af_add_cc_variable(tree, "err", ix=i_err)
42 call af_add_cc_variable(tree, "sol", ix=i_sol)
43 call af_add_cc_variable(tree, "tmp", ix=i_tmp)
44 call af_add_cc_variable(tree, "Dx", ix=i_gradx)
45 call af_add_cc_variable(tree, "eDx", ix=i_egradx)
46
47 ! Initialize tree
48 call af_init(tree, & ! Tree to initialize
49 box_size, & ! A box contains box_size**DIM cells
50 [1.25e-2_dp, 1.25e-2_dp], &
51 [box_size, box_size], &
52 coord=af_cyl) ! Cylindrical coordinates
53
54 call af_print_info(tree)
55
56 call system_clock(t_start, count_rate)
57 do
58 ! For each box, set the initial conditions
59 call af_loop_box(tree, set_init_cond)
60
61 ! This updates the refinement of the tree, by at most one level per call.
62 call af_adjust_refinement(tree, ref_routine, ref_info)
63
64 ! If no new boxes have been added, exit the loop
65 if (ref_info%n_add == 0) exit
66 end do
67 call system_clock(t_end, count_rate)
68
69 write(*,"(A,Es10.3,A)") " Wall-clock time generating AMR grid: ", &
70 (t_end-t_start) / real(count_rate,dp), " seconds"
71
72 call af_print_info(tree)
73
74 ! Set the multigrid options.
75 mg%i_phi = i_phi ! Solution variable
76 mg%i_rhs = i_rhs ! Right-hand side variable
77 mg%i_tmp = i_tmp ! Variable for temporary space
78 mg%sides_bc => sides_bc ! Method for boundary conditions Because we use
79
80 ! Initialize the multigrid options. This performs some basics checks and sets
81 ! default values where necessary.
82 ! This routine does not initialize the multigrid variables i_phi, i_rhs
83 ! and i_tmp. These variables will be initialized at the first call of mg_fas_fmg
84 call mg_init(tree, mg)
85
86 print *, "Multigrid iteration | max residual | max rel. error"
87 call system_clock(t_start, count_rate)
88
89 do mg_iter = 1, n_iterations
90 ! Perform a FAS-FMG cycle (full approximation scheme, full multigrid). The
91 ! third argument controls whether the residual is stored in i_tmp. The
92 ! fourth argument controls whether to improve the current solution.
93 call mg_fas_fmg(tree, mg, .true., mg_iter>1)
94
95 ! Compute the error compared to the analytic solution
96 call af_loop_box(tree, set_err)
97
98 ! Determine the minimum and maximum residual and error
99 call af_tree_min_cc(tree, i_tmp, residu(1))
100 call af_tree_max_cc(tree, i_tmp, residu(2))
101 call af_tree_min_cc(tree, i_err, anal_err(1))
102 call af_tree_max_cc(tree, i_err, anal_err(2))
103 call af_tree_max_cc(tree, i_phi, max_val)
104 anal_err = anal_err / max_val
105
106 write(*,"(I8,2Es14.5)") mg_iter, maxval(abs(residu)), &
107 maxval(abs(anal_err))
108
109 write(fname, "(A,I0)") "output/poisson_cyl_analytic_", mg_iter
110 call af_write_silo(tree, trim(fname))
111 end do
112 call system_clock(t_end, count_rate)
113
114 write(*, "(A,I0,A,E10.3,A)") &
115 " Wall-clock time after ", n_iterations, &
116 " iterations: ", (t_end-t_start) / real(count_rate, dp), &
117 " seconds"
118
119 ! This call is not really necessary here, but cleaning up the data in a tree
120 ! is important if your program continues with other tasks.
121 call af_destroy(tree)
122
123contains
124
125 ! Set refinement flags for box
126 subroutine ref_routine(box, cell_flags)
127 type(box_t), intent(in) :: box
128 integer, intent(out) :: cell_flags(box%n_cell, box%n_cell)
129 integer :: i, j, nc
130 real(dp) :: crv, dr2
131
132 nc = box%n_cell
133 dr2 = maxval(box%dr)**2
134
135 ! Compute the "curvature" in phi
136 do j = 1, nc
137 do i = 1, nc
138 crv = dr2 * abs(box%cc(i, j, i_rhs))
139
140 ! And refine if it exceeds a threshold
141 if (crv > 1e-1_dp .and. box%lvl < 10) then
142 cell_flags(i, j) = af_do_ref
143 else
144 cell_flags(i, j) = af_keep_ref
145 end if
146 end do
147 end do
148 end subroutine ref_routine
149
150 ! This routine sets the initial conditions for each box
151 subroutine set_init_cond(box)
152 type(box_t), intent(inout) :: box
153 integer :: i, j, nc
154 real(dp) :: rz(2), gradx
155
156 nc = box%n_cell
157
158 do j = 0, nc+1
159 do i = 0, nc+1
160 rz = af_r_cc(box, [i,j])
161
162 ! Gaussian source term
163 box%cc(i, j, i_rhs) = - q_source * &
164 exp(-sum((rz - rz_source)**2) / (2 * sigma**2)) / &
165 (sigma**3 * sqrt(2 * pi)**3 * epsilon_source)
166
167 box%cc(i, j, i_sol) = solution(rz)
168 end do
169 end do
170
171 do j = 1, nc
172 do i = 1, nc
173 gradx = 0.5_dp * (box%cc(i+1, j, i_sol) - &
174 box%cc(i-1, j, i_sol)) / box%dr(1)
175 box%cc(i, j, i_gradx) = gradx
176 end do
177 end do
178
179 end subroutine set_init_cond
180
181 real(dp) function solution(rz_in)
182 real(dp), intent(in) :: rz_in(2)
183 real(dp) :: d, erf_d
184
185 d = norm2(rz_in - rz_source)
186
187 ! Take limit when d -> 0
188 if (d < sqrt(epsilon(1.0d0))) then
189 erf_d = sqrt(2.0_dp/pi) / sigma
190 else
191 erf_d = erf(d * sqrt(0.5_dp) / sigma) / d
192 end if
193
194 solution = erf_d * q_source / (4 * pi * epsilon_source)
195 end function solution
196
197 ! Compute error compared to the analytic solution
198 subroutine set_err(box)
199 type(box_t), intent(inout) :: box
200 integer :: i, j, nc
201 real(dp) :: rz(2), gradx
202
203 nc = box%n_cell
204 do j = 1, nc
205 do i = 1, nc
206 rz = af_r_cc(box, [i,j])
207 box%cc(i, j, i_err) = box%cc(i, j, i_phi) - box%cc(i, j, i_sol)
208
209 gradx = 0.5_dp * (box%cc(i+1, j, i_phi) - &
210 box%cc(i-1, j, i_phi)) / box%dr(1)
211 box%cc(i, j, i_egradx) = gradx - box%cc(i, j, i_gradx)
212 end do
213 end do
214 end subroutine set_err
215
216 ! This routine sets boundary conditions for a box
217 subroutine sides_bc(box, nb, iv, coords, bc_val, bc_type)
218 type(box_t), intent(in) :: box
219 integer, intent(in) :: nb
220 integer, intent(in) :: iv
221 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
222 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
223 integer, intent(out) :: bc_type
224 integer :: n
225
226 if (nb == af_neighb_lowx) then
227 ! On the axis, apply Neumann zero conditions
228 bc_type = af_bc_neumann
229 bc_val = 0.0_dp
230 else
231 ! We use dirichlet boundary conditions
232 bc_type = af_bc_dirichlet
233
234 ! Below the solution is specified in the approriate ghost cells
235 do n = 1, box%n_cell**(ndim-1)
236 bc_val(n) = solution(coords(:, n))
237 end do
238 end if
239 end subroutine sides_bc
240
241end program poisson_cyl_analytic
Module which contains all Afivo modules, so that a user does not have to include them separately.
Definition: m_af_all.f90:3