Afivo 0.3
All Classes Namespaces Functions Variables Pages
dielectric_surface.f90

Example showing how to include a dielectric surface.

Example showing how to include a dielectric surface

1#include "../src/cpp_macros.h"
2!> \example dielectric_surface.f90
3!>
4!> Example showing how to include a dielectric surface
5program dielectric_surface
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_tmp
15 integer :: i_eps
16 integer :: i_fld_norm_fc
17 integer :: i_fld_cc(NDIM)
18 integer :: i_fld_fc
19 real(dp) :: fac = 1.0_dp
20
21 ! The dielectric constant used in this example
22 double precision, parameter :: epsilon_high = 2.0_dp
23
24 ! Where the interface is located
25 real(dp), parameter :: interface_location = 0.25_dp
26 ! Along which dimension the interface occurs
27 integer, parameter :: interface_dimension = 1
28
29 type(af_t) :: tree
30 type(ref_info_t) :: ref_info
31 type(surfaces_t) :: dielectric
32 integer :: n, mg_iter
33 real(dp) :: residu
34 character(len=100) :: fname
35 type(mg_t) :: mg
36 integer, allocatable :: ref_links(:, :)
37
38 call af_add_cc_variable(tree, "phi", ix=i_phi)
39 call af_add_cc_variable(tree, "rhs", ix=i_rhs)
40 call af_add_cc_variable(tree, "tmp", ix=i_tmp)
41 call af_add_cc_variable(tree, "eps", ix=i_eps)
42
43 ! Include both cell-centered and face-centered electric field for testing
44 call af_add_cc_variable(tree, "fld_cc_x", ix=i_fld_cc(1))
45#if NDIM > 1
46 call af_add_cc_variable(tree, "fld_cc_y", ix=i_fld_cc(2))
47#endif
48#if NDIM == 3
49 call af_add_cc_variable(tree, "fld_cc_z", ix=i_fld_cc(3))
50#endif
51 call af_add_cc_variable(tree, "fld_norm_fc", ix=i_fld_norm_fc)
52 call af_add_fc_variable(tree, "fld_fc", i_fld_fc)
53
54 call af_set_cc_methods(tree, i_eps, af_bc_neumann_zero, &
55 af_gc_prolong_copy, prolong=af_prolong_zeroth)
56 call af_set_cc_methods(tree, i_rhs, af_bc_neumann_zero)
57
58 ! Initialize tree
59 call af_init(tree, & ! Tree to initialize
60 box_size, & ! A box contains box_size**DIM cells
61 [dtimes(1.0_dp)], &
62 [dtimes(box_size)])
63
64 do n = 1, 2
65 call af_adjust_refinement(tree, ref_routine, ref_info)
66 if (ref_info%n_add == 0) exit
67 end do
68
69 call af_loop_box(tree, set_init_cond)
70
71 call surface_initialize(tree, i_eps, dielectric, 1)
72
73 call surface_set_values(tree, dielectric, 1, sigma_function)
74
75 do n = 1, 3
76 call surface_get_refinement_links(dielectric, ref_links)
77 call af_adjust_refinement(tree, ref_random, ref_info, ref_links=ref_links)
78 call surface_update_after_refinement(tree, dielectric, ref_info)
79 end do
80
81 call surface_surface_charge_to_rhs(tree, dielectric, 1, i_rhs, fac)
82
83 tree%mg_i_eps = i_eps
84 mg%i_phi = i_phi
85 mg%i_rhs = i_rhs
86 mg%i_tmp = i_tmp
87 mg%sides_bc => bc_phi
88
89 call mg_init(tree, mg)
90
91 do mg_iter = 1, n_iterations
92 call mg_fas_fmg(tree, mg, .true., mg_iter>1)
93 call af_loop_box(tree, compute_fields)
94 call surface_correct_field_fc(tree, dielectric, 1, i_fld_fc, i_phi, -fac)
95 call af_loop_box(tree, compute_field_norms)
96
97 ! Determine the minimum and maximum residual and error
98 call af_tree_maxabs_cc(tree, i_tmp, residu)
99 write(*, "(I8,Es14.5)") mg_iter, residu
100
101 write(fname, "(A,I0)") "output/dielectric_surface_" // &
102 dimname // "_", mg_iter
103 call af_write_silo(tree, trim(fname))
104 end do
105
106contains
107
108 ! Set refinement flags for box
109 subroutine ref_routine(box, cell_flags)
110 type(box_t), intent(in) :: box
111 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
112
113 if (box%lvl < 3) then
114 cell_flags(dtimes(:)) = af_do_ref
115 else
116 cell_flags(dtimes(:)) = af_keep_ref
117 end if
118 end subroutine ref_routine
119
120 subroutine ref_random(box, cell_flags)
121 type(box_t), intent(in) :: box ! A list of all boxes in the tree
122 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
123 real(dp) :: rr
124
125 ! Draw a [0, 1) random number
126 call random_number(rr)
127
128 if (rr < 0.5_dp**ndim .and. box%lvl < 6) then
129 cell_flags = af_do_ref
130 else if (box%lvl > 3) then
131 cell_flags = af_rm_ref
132 else
133 cell_flags = af_keep_ref
134 end if
135 end subroutine ref_random
136
137 ! This routine sets the initial conditions for each box
138 subroutine set_init_cond(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(0,nc+1)
146 rr = af_r_cc(box, [ijk])
147
148 ! Change epsilon in part of the domain
149 if (rr(interface_dimension) < interface_location) then
150 box%cc(ijk, i_eps) = epsilon_high
151 else
152 box%cc(ijk, i_eps) = 1.0_dp
153 end if
154
155 ! box%cc(IJK, i_rhs) = 10 * epsilon_high * &
156 ! exp(-100 * sum((rr - 0.75_dp)**2))
157 end do; close_do
158
159 end subroutine set_init_cond
160
161 ! This routine sets boundary conditions for a box
162 subroutine bc_phi(box, nb, iv, coords, bc_val, bc_type)
163 type(box_t), intent(in) :: box
164 integer, intent(in) :: nb
165 integer, intent(in) :: iv
166 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
167 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
168 integer, intent(out) :: bc_type
169 integer :: nc
170
171 nc = box%n_cell
172
173 ! Below the solution is specified in the approriate ghost cells
174 select case (nb)
175 case (2 * interface_dimension - 1)
176 bc_type = af_bc_dirichlet
177 bc_val = 0.0_dp
178 case (2 * interface_dimension)
179 bc_type = af_bc_dirichlet
180 bc_val = 1.0_dp
181 case default
182 bc_type = af_bc_neumann
183 bc_val = 0.0_dp
184 end select
185 end subroutine bc_phi
186
187 real(dp) function sigma_function(r)
188 real(dp), intent(in) :: r(NDIM)
189
190 ! Screen electric field outside dielectric
191 ! sigma_function = -epsilon_high/interface_location
192
193 ! Screen electric field inside dielectric
194 ! sigma_function = 1/(1 - interface_location)
195
196 ! Equal electric field on left and right
197 ! sigma_function = 1 - epsilon_high
198
199 ! No surface charge
200 sigma_function = 0.0_dp
201 end function sigma_function
202
203 subroutine compute_fields(box)
204 type(box_t), intent(inout) :: box
205 integer :: nc
206 real(dp) :: inv_dr(NDIM)
207
208 nc = box%n_cell
209 inv_dr = 1 / box%dr
210
211#if NDIM == 2
212 box%fc(1:nc+1, 1:nc, 1, i_fld_fc) = inv_dr(1) * &
213 (box%cc(0:nc, 1:nc, i_phi) - box%cc(1:nc+1, 1:nc, i_phi))
214 box%fc(1:nc, 1:nc+1, 2, i_fld_fc) = inv_dr(2) * &
215 (box%cc(1:nc, 0:nc, i_phi) - box%cc(1:nc, 1:nc+1, i_phi))
216 box%cc(1:nc, 1:nc, i_fld_cc(1)) = 0.5_dp * inv_dr(1) * &
217 (box%cc(0:nc-1, 1:nc, i_phi) - box%cc(2:nc+1, 1:nc, i_phi))
218 box%cc(1:nc, 1:nc, i_fld_cc(2)) = 0.5_dp * inv_dr(2) * &
219 (box%cc(1:nc, 0:nc-1, i_phi) - box%cc(1:nc, 2:nc+1, i_phi))
220#elif NDIM == 3
221 box%fc(1:nc+1, 1:nc, 1:nc, 1, i_fld_fc) = inv_dr(1) * &
222 (box%cc(0:nc, 1:nc, 1:nc, i_phi) - &
223 box%cc(1:nc+1, 1:nc, 1:nc, i_phi))
224 box%fc(1:nc, 1:nc+1, 1:nc, 2, i_fld_fc) = inv_dr(2) * &
225 (box%cc(1:nc, 0:nc, 1:nc, i_phi) - &
226 box%cc(1:nc, 1:nc+1, 1:nc, i_phi))
227 box%fc(1:nc, 1:nc, 1:nc+1, 3, i_fld_fc) = inv_dr(3) * &
228 (box%cc(1:nc, 1:nc, 0:nc, i_phi) - &
229 box%cc(1:nc, 1:nc, 1:nc+1, i_phi))
230#endif
231
232 end subroutine compute_fields
233
234 subroutine compute_field_norms(box)
235 type(box_t), intent(inout) :: box
236 integer :: nc
237
238 nc = box%n_cell
239
240#if NDIM == 2
241 box%cc(1:nc, 1:nc, i_fld_norm_fc) = 0.5_dp * sqrt(&
242 (box%fc(1:nc, 1:nc, 1, i_fld_fc) + &
243 box%fc(2:nc+1, 1:nc, 1, i_fld_fc))**2 + &
244 (box%fc(1:nc, 1:nc, 2, i_fld_fc) + &
245 box%fc(1:nc, 2:nc+1, 2, i_fld_fc))**2)
246#elif NDIM == 3
247 box%cc(1:nc, 1:nc, 1:nc, i_fld_norm_fc) = 0.5_dp * sqrt(&
248 (box%fc(1:nc, 1:nc, 1:nc, 1, i_fld_fc) + &
249 box%fc(2:nc+1, 1:nc, 1:nc, 1, i_fld_fc))**2 + &
250 (box%fc(1:nc, 1:nc, 1:nc, 2, i_fld_fc) + &
251 box%fc(1:nc, 2:nc+1, 1:nc, 2, i_fld_fc))**2 + &
252 (box%fc(1:nc, 1:nc, 1:nc, 3, i_fld_fc) + &
253 box%fc(1:nc, 1:nc, 2:nc+1, 3, i_fld_fc))**2)
254#endif
255 end subroutine compute_field_norms
256
257end program dielectric_surface
Module which contains all Afivo modules, so that a user does not have to include them separately.
Definition: m_af_all.f90:3