Afivo 0.3
All Classes Namespaces Functions Variables Pages
electrode_dielectric.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 electrode_dielectric.f90
3!>
4!> Example showing how to include a dielectric surface
5program electrode_dielectric
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_lsf
16 integer :: i_eps
17 integer :: i_field
18 integer :: i_field_norm
19
20 type(af_t) :: tree
21 type(ref_info_t) :: ref_info
22 integer :: n, mg_iter, coord, n_args
23 real(dp) :: residu
24 character(len=100) :: fname
25 type(mg_t) :: mg
26
27 if (ndim /= 2) error stop "Example only set up for 2D"
28
29 call af_add_cc_variable(tree, "phi", ix=i_phi)
30 call af_add_cc_variable(tree, "rhs", ix=i_rhs)
31 call af_add_cc_variable(tree, "tmp", ix=i_tmp)
32 call af_add_cc_variable(tree, "lsf", ix=i_lsf)
33 call af_add_cc_variable(tree, "eps", ix=i_eps)
34 call af_add_cc_variable(tree, "field_norm", ix=i_field_norm)
35 call af_add_fc_variable(tree, "field", ix=i_field)
36
37 call af_set_cc_methods(tree, i_lsf, funcval=set_lsf)
38 call af_set_cc_methods(tree, i_rhs, af_bc_neumann_zero)
39
40 ! If an argument is given, switch to cylindrical coordinates in 2D
41 n_args = command_argument_count()
42 if (ndim == 2 .and. n_args == 1) then
43 coord = af_cyl
44 else
45 coord = af_xyz
46 end if
47
48 ! Initialize tree
49 call af_init(tree, & ! Tree to initialize
50 box_size, & ! A box contains box_size**DIM cells
51 [dtimes(1.0_dp)], &
52 4 * [dtimes(box_size)], &
53 coord=coord)
54
55 do n = 1, 4
56 call af_adjust_refinement(tree, ref_routine, ref_info)
57 if (ref_info%n_add == 0) exit
58 end do
59
60 tree%mg_i_lsf = i_lsf
61 tree%mg_i_eps = i_eps
62 mg%i_phi = i_phi
63 mg%i_rhs = i_rhs
64 mg%i_tmp = i_tmp
65 mg%sides_bc => af_bc_dirichlet_zero
66 mg%lsf => get_lsf
67 mg%lsf_boundary_value = 1.0_dp
68
69 call mg_init(tree, mg)
70
71 do mg_iter = 1, n_iterations
72 call mg_fas_fmg(tree, mg, .true., mg_iter>1)
73 call mg_compute_phi_gradient(tree, mg, i_field, 1.0_dp, i_field_norm)
74
75 ! Determine the minimum and maximum residual and error
76 call af_tree_maxabs_cc(tree, i_tmp, residu)
77 write(*, "(I8,Es14.5)") mg_iter, residu
78
79 write(fname, "(A,I0)") "output/electrode_dielectric_" // &
80 dimname // "_", mg_iter
81 call af_write_silo(tree, trim(fname))
82 end do
83
84contains
85
86 ! Set refinement flags for box
87 subroutine ref_routine(box, cell_flags)
88 type(box_t), intent(in) :: box
89 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
90
91 if (box%lvl < 9 - 2 * ndim .and. box%r_min(2) < 0.5_dp) then
92 cell_flags(dtimes(:)) = af_do_ref
93 else
94 cell_flags(dtimes(:)) = af_keep_ref
95 end if
96 end subroutine ref_routine
97
98 ! This routine sets the level set function
99 subroutine set_lsf(box, iv)
100 type(box_t), intent(inout) :: box
101 integer, intent(in) :: iv
102 real(dp) :: rr(NDIM)
103 integer :: IJK, nc
104
105 nc = box%n_cell
106
107 do kji_do(0,nc+1)
108 rr = af_r_cc(box, [ijk])
109 box%cc(ijk, iv) = get_lsf(rr)
110 if (rr(1) < 0.5_dp) then
111 box%cc(ijk, i_eps) = 20.0_dp
112 else
113 box%cc(ijk, i_eps) = 1.0_dp
114 end if
115 end do; close_do
116 end subroutine set_lsf
117
118 real(dp) function get_lsf(r)
119 real(dp), intent(in) :: r(NDIM)
120 real(dp) :: r0(NDIM), r1(NDIM)
121 real(dp) :: radius
122
123 ! Start and end point of line
124 r0(:) = 0.3_dp
125 r1(:) = 0.5_dp
126 radius = 0.02_dp
127 get_lsf = dist_vec_line(r, r0, r1, ndim) - radius
128 end function get_lsf
129
130 !> Compute distance vector between point and its projection onto a line
131 !> between r0 and r1
132 function dist_vec_line(r, r0, r1, n_dim) result(dist)
133 integer, intent(in) :: n_dim
134 real(dp), intent(in) :: r(n_dim), r0(n_dim), r1(n_dim)
135 real(dp) :: dist_vec(n_dim), frac
136 real(dp) :: dist, line_len2
137
138 line_len2 = sum((r1 - r0)**2)
139 frac = sum((r - r0) * (r1 - r0))
140
141 if (frac <= 0.0_dp) then
142 dist_vec = r - r0
143 else if (frac >= line_len2) then
144 dist_vec = r - r1
145 else
146 dist_vec = r - (r0 + frac/line_len2 * (r1 - r0))
147 end if
148
149 dist = norm2(dist_vec)
150 end function dist_vec_line
151
152end program electrode_dielectric
Module which contains all Afivo modules, so that a user does not have to include them separately.
Definition: m_af_all.f90:3