Afivo 0.3
All Classes Namespaces Functions Variables Pages
poisson_cyl_dielectric.f90

Example showing how to use m_af_multigrid in cylindrical coordinates with an abrubt change in "eps", and compare with an analytic solution.

Example showing how to use m_af_multigrid in cylindrical coordinates with an abrubt change in "eps", and compare with an analytic solution.

1!> \example poisson_cyl_dielectric.f90
2!>
3!> Example showing how to use m_af_multigrid in cylindrical coordinates with an abrubt
4!> change in "eps", and compare with an analytic solution.
5program poisson_cyl_dielectric
6 use m_af_all
8
9 implicit none
10
11 integer, parameter :: box_size = 8
12 integer, parameter :: n_iterations = 10
13 integer :: i_phi
14 integer :: i_rhs
15 integer :: i_err
16 integer :: i_tmp
17 integer :: i_eps
18
19 type(af_t) :: tree
20 type(ref_info_t) :: ref_info
21 integer :: mg_iter
22 real(dp) :: residu(2), anal_err(2)
23 character(len=100) :: fname
24 type(mg_t) :: mg
25 type(gauss_t) :: gs
26 integer :: count_rate,t_start, t_end
27
28 print *, "Running poisson_cyl_dielectric"
29 print *, "Number of threads", af_get_max_threads()
30
31 ! The manufactured solution exists of two Gaussians, which are stored in gs
32 call gauss_init(gs, [1.0_dp], [0.1_dp], &
33 reshape([0.0_dp, 0.25_dp], [2,1]))
34
35 call af_add_cc_variable(tree, "phi", ix=i_phi)
36 call af_add_cc_variable(tree, "rhs", ix=i_rhs)
37 call af_add_cc_variable(tree, "err", ix=i_err)
38 call af_add_cc_variable(tree, "tmp", ix=i_tmp)
39 call af_add_cc_variable(tree, "eps", ix=i_eps)
40
41 ! Initialize tree
42 call af_init(tree, & ! Tree to initialize
43 box_size, & ! A box contains box_size**DIM cells
44 [1.0_dp, 1.0_dp], &
45 [box_size, box_size], &
46 coord=af_cyl) ! Cylindrical coordinates
47
48 call af_print_info(tree)
49
50 call system_clock(t_start, count_rate)
51 do
52 ! For each box, set the initial conditions
53 call af_loop_box(tree, set_init_cond)
54
55 ! This updates the refinement of the tree, by at most one level per call.
56 call af_adjust_refinement(tree, ref_routine, ref_info)
57
58 ! If no new boxes have been added, exit the loop
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 ! Set the multigrid options.
69 tree%mg_i_eps = i_eps
70 mg%i_phi = i_phi ! Solution variable
71 mg%i_rhs = i_rhs ! Right-hand side variable
72 mg%i_tmp = i_tmp ! Variable for temporary space
73 mg%sides_bc => sides_bc ! Method for boundary conditions Because we use
74
75 ! Initialize the multigrid options. This performs some basics checks and sets
76 ! default values where necessary.
77 ! This routine does not initialize the multigrid variables i_phi, i_rhs
78 ! and i_tmp. These variables will be initialized at the first call of mg_fas_fmg
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 ! Perform a FAS-FMG cycle (full approximation scheme, full multigrid). The
86 ! third argument controls whether the residual is stored in i_tmp. The
87 ! fourth argument controls whether to improve the current solution.
88 call mg_fas_fmg(tree, mg, .true., mg_iter>1)
89
90 ! Compute the error compared to the analytic solution
91 call af_loop_box(tree, set_err)
92
93 ! Determine the minimum and maximum residual and error
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/poisson_cyl_dielectric_", mg_iter
102 call af_write_silo(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 ! This call is not really necessary here, but cleaning up the data in a tree
112 ! is important if your program continues with other tasks.
113 call af_destroy(tree)
114
115contains
116
117 ! Set refinement flags for box
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 ! Compute the "curvature" in phi
128 do j = 1, nc
129 do i = 1, nc
130 crv = dr2 * abs(box%cc(i, j, i_rhs)) / box%cc(i, j, i_eps)
131
132 ! And refine if it exceeds a threshold
133 if (crv > 1.0e-3_dp) then
134 cell_flags(i, j) = af_do_ref
135 else
136 cell_flags(i, j) = af_keep_ref
137 end if
138 end do
139 end do
140 end subroutine ref_routine
141
142 ! This routine sets the initial conditions for each box
143 subroutine set_init_cond(box)
144 type(box_t), intent(inout) :: box
145 integer :: i, j, nc
146 real(dp) :: rz(2), grad(2), qbnd, tmp
147
148 nc = box%n_cell
149 box%cc(:, :, i_phi) = 0
150
151 do j = 0, nc+1
152 do i = 0, nc+1
153 rz = af_r_cc(box, [i,j])
154
155 ! Change epsilon in part of the domain
156 if (rz(1) < 0.5_dp .and. rz(2) < 0.5_dp) then
157 box%cc(i, j, i_eps) = 100.0_dp
158 else
159 box%cc(i, j, i_eps) = 1.0_dp
160 end if
161
162 ! Partially compute the right-hand side (see below)
163 box%cc(i, j, i_rhs) = gauss_laplacian_cyl(gs, rz) * box%cc(i, j, i_eps)
164 end do
165 end do
166
167 ! We have to place surface charges where epsilon has a jump, this is first
168 ! done in the r-direction
169 do j = 1, nc
170 do i = 0, nc
171 rz = af_rr_cc(box, [i + 0.5_dp, real(j, dp)])
172
173 ! Determine amount of charge
174 call gauss_gradient(gs, rz, grad)
175 qbnd = (box%cc(i+1, j, i_eps) - box%cc(i, j, i_eps)) * &
176 grad(1) / box%dr(1)
177
178 ! Place surface charge weighted with eps
179 tmp = box%cc(i+1, j, i_eps) / &
180 (box%cc(i, j, i_eps) + box%cc(i+1, j, i_eps))
181 box%cc(i+1, j, i_rhs) = box%cc(i+1, j, i_rhs) + tmp * qbnd
182 box%cc(i, j, i_rhs) = box%cc(i, j, i_rhs) + (1-tmp) * qbnd
183 end do
184 end do
185
186 ! Set surface charge in z-direction
187 do j = 0, nc
188 do i = 1, nc
189 rz = af_rr_cc(box, [real(i, dp), j + 0.5_dp])
190
191 ! Determine amount of charge
192 call gauss_gradient(gs, rz, grad)
193 qbnd = (box%cc(i, j+1, i_eps) - box%cc(i, j, i_eps)) * &
194 grad(2) / box%dr(2)
195
196 ! Place surface charge weighted with eps
197 tmp = box%cc(i, j+1, i_eps) / &
198 (box%cc(i, j, i_eps) + box%cc(i, j+1, i_eps))
199 box%cc(i, j+1, i_rhs) = box%cc(i, j+1, i_rhs) + tmp * qbnd
200 box%cc(i, j, i_rhs) = box%cc(i, j, i_rhs) + (1-tmp) * qbnd
201 end do
202 end do
203 end subroutine set_init_cond
204
205 ! Compute error compared to the analytic solution
206 subroutine set_err(box)
207 type(box_t), intent(inout) :: box
208 integer :: i, j, nc
209 real(dp) :: rz(2)
210
211 nc = box%n_cell
212 do j = 1, nc
213 do i = 1, nc
214 rz = af_r_cc(box, [i,j])
215 box%cc(i, j, i_err) = box%cc(i, j, i_phi) - gauss_value(gs, rz)
216 end do
217 end do
218 end subroutine set_err
219
220 ! This routine sets boundary conditions for a box
221 subroutine sides_bc(box, nb, iv, coords, bc_val, bc_type)
222 type(box_t), intent(in) :: box
223 integer, intent(in) :: nb
224 integer, intent(in) :: iv
225 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
226 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
227 integer, intent(out) :: bc_type
228 integer :: n
229
230 if (nb == af_neighb_lowx) then
231 ! On the axis, apply Neumann zero conditions
232 bc_type = af_bc_neumann
233 bc_val = 0.0_dp
234 else
235 ! We use dirichlet boundary conditions
236 bc_type = af_bc_dirichlet
237
238 ! Below the solution is specified in the approriate ghost cells
239 do n = 1, box%n_cell**(ndim-1)
240 bc_val(n) = gauss_value(gs, coords(:, n))
241 end do
242 end if
243 end subroutine sides_bc
244
245end program poisson_cyl_dielectric
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