Afivo 0.3
All Classes Namespaces Functions Variables Pages
two_electrodes.f90

Test with two electrodes.

Test with two electrodes

1#include "../src/cpp_macros.h"
2!> \example two_electrodes.f90
3!>
4!> Test with two electrodes
5program two_electrodes
6 use m_af_all
7 use m_config
8 use omp_lib
9
10 implicit none
11
12 integer :: box_size = 8
13 integer :: n_iterations = 10
14 integer :: max_refine_level = 7
15 integer :: min_refine_level = 5
16 integer :: i_phi
17 integer :: i_rhs
18 integer :: i_tmp
19 integer :: i_lsf
20 integer :: i_field
21 integer :: i_field_norm
22
23 type(af_t) :: tree
24 type(ref_info_t) :: ref_info
25 integer :: n, mg_iter
26 real(dp) :: residu
27 character(len=150) :: fname
28 type(mg_t) :: mg
29 real(dp) :: mem_limit_gb = 8.0_dp
30 real(dp) :: t0, t_sum
31
32 real(dp) :: rod_r0(NDIM), rod_r1(NDIM), rod_radius
33 real(dp) :: sphere_r0(NDIM), sphere_radius
34
35 rod_r0 = 0.5_dp
36 rod_r0(ndim) = 0.7_dp
37 rod_r1 = 0.5_dp
38 rod_r1(ndim) = 1.0_dp
39 rod_radius = 0.05_dp
40
41 sphere_r0 = 0.5_dp
42 sphere_r0(ndim) = 0.0_dp
43 sphere_radius = 0.25_dp
44
45 call af_add_cc_variable(tree, "phi", ix=i_phi)
46 call af_add_cc_variable(tree, "rhs", ix=i_rhs)
47 call af_add_cc_variable(tree, "tmp", ix=i_tmp)
48 call af_add_cc_variable(tree, "lsf", ix=i_lsf)
49 call af_add_cc_variable(tree, "field_norm", ix=i_field_norm)
50 call af_add_fc_variable(tree, "field", ix=i_field)
51
52 call af_set_cc_methods(tree, i_lsf, funcval=set_lsf)
53
54 tree%mg_i_lsf = i_lsf
55 mg%i_phi = i_phi
56 mg%i_rhs = i_rhs
57 mg%i_tmp = i_tmp
58 mg%sides_bc => af_bc_dirichlet_zero
59 mg%lsf_boundary_value = 1.0_dp
60 mg%lsf => get_lsf
61 mg%lsf_length_scale = 1e-2_dp
62 mg%lsf_boundary_function => get_boundary_value
63 mg%sides_bc => bc_sides
64
65 ! Initialize tree
66 call af_init(tree, & ! Tree to initialize
67 box_size, & ! A box contains box_size**DIM cells
68 [dtimes(1.0_dp)], &
69 [dtimes(box_size)], &
70 mem_limit_gb=mem_limit_gb)
71
72 call af_refine_up_to_lvl(tree, min_refine_level)
73 call mg_init(tree, mg)
74
75 do n = 1, 100
76 call af_adjust_refinement(tree, ref_routine, ref_info)
77 if (ref_info%n_add == 0) exit
78 end do
79
80 call af_print_info(tree)
81
82 t_sum = 0
83 write(*, "(A,A4,4A14)") "# ", "iter", "residu"
84
85 do mg_iter = 1, n_iterations
86 t0 = omp_get_wtime()
87 call mg_fas_fmg(tree, mg, .true., mg_iter>1)
88 t_sum = t_sum + omp_get_wtime() - t0
89
90 call mg_compute_phi_gradient(tree, mg, i_field, 1.0_dp, i_field_norm)
91 call af_loop_box(tree, set_field_to_zero_inside)
92
93 ! Determine the minimum and maximum residual and error
94 call af_tree_maxabs_cc(tree, i_tmp, residu)
95 write(*, "(A,i4,4E14.5)") "# ", mg_iter, residu
96
97 write(fname, "(A,I0)") "output/two_electrodes_" // dimname // "_", mg_iter
98 call af_write_silo(tree, trim(fname))
99 end do
100
101 write(*, "(A,E14.5)") " seconds_per_FMG_cycle", t_sum/n_iterations
102
103 call af_stencil_print_info(tree)
104
105contains
106
107 ! Set refinement flags for box
108 subroutine ref_routine(box, cell_flags)
109 type(box_t), intent(in) :: box
110 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
111 integer :: IJK, nc
112 real(dp) :: r(NDIM), r_ref(NDIM)
113
114 nc = box%n_cell
115 cell_flags(dtimes(:)) = af_keep_ref
116
117 ! Refine around this point
118 r_ref = rod_r0
119 r_ref(ndim) = r_ref(ndim) - rod_radius
120
121 if (box%lvl < max_refine_level) then
122 do kji_do(1, nc)
123 r = af_r_cc(box, [ijk])
124 if (norm2(r - r_ref) < 2 * rod_radius) then
125 cell_flags(dtimes(:)) = af_do_ref
126 end if
127 end do; close_do
128 end if
129 end subroutine ref_routine
130
131 ! This routine sets the level set function
132 subroutine set_lsf(box, iv)
133 type(box_t), intent(inout) :: box
134 integer, intent(in) :: iv
135 integer :: IJK, nc
136 real(dp) :: rr(NDIM)
137
138 nc = box%n_cell
139
140 do kji_do(0,nc+1)
141 rr = af_r_cc(box, [ijk])
142 box%cc(ijk, iv) = get_lsf(rr)
143 end do; close_do
144 end subroutine set_lsf
145
146 ! The value at the lsf boundaries
147 real(dp) function get_boundary_value(r)
148 real(dp), intent(in) :: r(NDIM)
149
150 if (r(ndim) > 0.5_dp) then
151 get_boundary_value = 1.0_dp
152 else
153 get_boundary_value = 0.0_dp
154 end if
155 end function get_boundary_value
156
157 real(dp) function get_lsf(r)
158 real(dp), intent(in) :: r(NDIM)
159
160 if (r(ndim) > 0.5_dp) then
161 get_lsf = gm_dist_line(r, rod_r0, rod_r1, ndim) - rod_radius
162 else
163 get_lsf = norm2(r - sphere_r0) - sphere_radius
164 end if
165 end function get_lsf
166
167 !> To set boundary conditions at the sides of the domain
168 subroutine bc_sides(box, nb, iv, coords, bc_val, bc_type)
169 type(box_t), intent(in) :: box
170 integer, intent(in) :: nb
171 integer, intent(in) :: iv
172 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
173 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
174 integer, intent(out) :: bc_type
175
176 bc_type = af_bc_dirichlet
177
178 if (af_neighb_dim(nb) == ndim) then
179 if (af_neighb_low(nb)) then
180 bc_type = af_bc_dirichlet
181 bc_val = 0.0_dp
182 else
183 bc_type = af_bc_dirichlet
184 bc_val = 1.0_dp
185 end if
186 else
187 bc_type = af_bc_neumann
188 bc_val = 0.0_dp
189 end if
190 end subroutine bc_sides
191
192 !> To set the electric field to zero inside the object
193 subroutine set_field_to_zero_inside(box)
194 type(box_t), intent(inout) :: box
195
196 where (box%cc(dtimes(:), i_lsf) < 0)
197 box%cc(dtimes(:), i_field_norm) = 0
198 end where
199 end subroutine set_field_to_zero_inside
200
201 !> Compute distance vector between point and its projection onto a line
202 !> between r0 and r1
203 subroutine gm_dist_vec_line(r, r0, r1, n_dim, dist_vec, frac)
204 integer, intent(in) :: n_dim
205 real(dp), intent(in) :: r(n_dim), r0(n_dim), r1(n_dim)
206 real(dp), intent(out) :: dist_vec(n_dim)
207 real(dp), intent(out) :: frac !< Fraction [0,1] along line
208 real(dp) :: line_len2
209
210 line_len2 = sum((r1 - r0)**2)
211 frac = sum((r - r0) * (r1 - r0))
212
213 if (frac <= 0.0_dp) then
214 frac = 0.0_dp
215 dist_vec = r - r0
216 else if (frac >= line_len2) then
217 frac = 1.0_dp
218 dist_vec = r - r1
219 else
220 dist_vec = r - (r0 + frac/line_len2 * (r1 - r0))
221 frac = sqrt(frac / line_len2)
222 end if
223 end subroutine gm_dist_vec_line
224
225 function gm_dist_line(r, r0, r1, n_dim) result(dist)
226 integer, intent(in) :: n_dim
227 real(dp), intent(in) :: r(n_dim), r0(n_dim), r1(n_dim)
228 real(dp) :: dist, dist_vec(n_dim), frac
229 call gm_dist_vec_line(r, r0, r1, n_dim, dist_vec, frac)
230 dist = norm2(dist_vec)
231 end function gm_dist_line
232
233end program
Module which contains all Afivo modules, so that a user does not have to include them separately.
Definition: m_af_all.f90:3
Module that allows working with a configuration file.
Definition: m_config.f90:5