Afivo 0.3
All Classes Namespaces Functions Variables Pages
poisson_dielectric.f90

Example showing how to include a dielectric object. Warning: the functionality is not fully ready.

Example showing how to include a dielectric object. Warning: the functionality is not fully ready

1#include "../src/cpp_macros.h"
2!> \example poisson_dielectric.f90
3!>
4!> Example showing how to include a dielectric object. Warning: the
5!> functionality is not fully ready
6program dielectric_test
7 use m_af_all
9
10 implicit none
11
12 integer, parameter :: box_size = 8
13 integer, parameter :: n_iterations = 10
14 integer :: i_phi
15 integer :: i_rhs
16 integer :: i_tmp
17 integer :: i_eps
18
19 ! The dielectric constant used in this example
20 double precision, parameter :: epsilon_high = 10.0_dp
21
22 type(af_t) :: tree
23 type(ref_info_t) :: ref_info
24 integer :: mg_iter
25 real(dp) :: residu(2)
26 character(len=100) :: fname
27 type(mg_t) :: mg
28 integer :: count_rate, t_start, t_end
29
30 print *, "****************************************"
31 print *, "Warning: functionality demonstrated here is not fully ready"
32 print *, "For large epsilon, convergence will probably be slow"
33 print *, "****************************************"
34 print *, "Number of threads", af_get_max_threads()
35
36 call af_add_cc_variable(tree, "phi", ix=i_phi)
37 call af_add_cc_variable(tree, "rhs", ix=i_rhs)
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 [dtimes(1.0_dp)], &
45 [dtimes(box_size)])
46
47 call af_print_info(tree)
48
49 call system_clock(t_start, count_rate)
50 do
51 ! For each box, set the initial conditions
52 call af_loop_box(tree, set_init_cond)
53
54 ! This updates the refinement of the tree, by at most one level per call.
55 call af_adjust_refinement(tree, ref_routine, ref_info)
56
57 ! If no new boxes have been added, exit the loop
58 if (ref_info%n_add == 0) exit
59 end do
60 call system_clock(t_end, count_rate)
61
62 ! Average epsilon on coarse grids. In the future, it could be better to define
63 ! epsilon on cell faces, and to perform this restriction in a matrix fashion:
64 ! A_coarse = M_restrict * A_fine * M_prolong (A = matrix operator, M = matrix)
65 call af_restrict_tree(tree, [i_eps])
66
67 write(*,"(A,Es10.3,A)") " Wall-clock time generating AMR grid: ", &
68 (t_end-t_start) / real(count_rate,dp), " seconds"
69
70 call af_print_info(tree)
71
72 ! Set the multigrid options.
73 tree%mg_i_eps = i_eps
74 mg%i_phi = i_phi ! Solution variable
75 mg%i_rhs = i_rhs ! Right-hand side variable
76 mg%i_tmp = i_tmp ! Variable for temporary space
77 mg%sides_bc => sides_bc
78
79 ! Initialize the multigrid options. This performs some basics checks and sets
80 ! default values where necessary.
81 ! This routine does not initialize the multigrid variables i_phi, i_rhs
82 ! and i_tmp. These variables will be initialized at the first call of mg_fas_fmg
83 call mg_init(tree, mg)
84
85 print *, "Multigrid iteration | max residual | max error"
86 call system_clock(t_start, count_rate)
87
88 do mg_iter = 1, n_iterations
89 ! Perform a FAS-FMG cycle (full approximation scheme, full multigrid). The
90 ! third argument controls whether the residual is stored in i_tmp. The
91 ! fourth argument controls whether to improve the current solution.
92 call mg_fas_fmg(tree, mg, .true., mg_iter>1)
93
94 ! Determine the minimum and maximum residual and error
95 call af_tree_min_cc(tree, i_tmp, residu(1))
96 call af_tree_max_cc(tree, i_tmp, residu(2))
97 write(*,"(I8,Es14.5)") mg_iter, maxval(abs(residu))
98
99 write(fname, "(A,I0)") "output/dielectric_" // dimname // "_", mg_iter
100 call af_write_silo(tree, trim(fname))
101 end do
102 call system_clock(t_end, count_rate)
103
104 write(*, "(A,I0,A,E10.3,A)") &
105 " Wall-clock time after ", n_iterations, &
106 " iterations: ", (t_end-t_start) / real(count_rate, dp), &
107 " seconds"
108
109 ! This call is not really necessary here, but cleaning up the data in a tree
110 ! is important if your program continues with other tasks.
111 call af_destroy(tree)
112
113contains
114
115 ! Set refinement flags for box
116 subroutine ref_routine(box, cell_flags)
117 type(box_t), intent(in) :: box
118 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
119 real(dp) :: eps_min, eps_max
120
121 eps_min = minval(box%cc(dtimes(:), i_eps))
122 eps_max = maxval(box%cc(dtimes(:), i_eps))
123
124 if ((box%lvl < 5 .and. eps_max > eps_min) .or. box%lvl < 2) then
125 cell_flags(dtimes(:)) = af_do_ref
126 else
127 cell_flags(dtimes(:)) = af_keep_ref
128 end if
129 end subroutine ref_routine
130
131 ! This routine sets the initial conditions for each box
132 subroutine set_init_cond(box)
133 type(box_t), intent(inout) :: box
134 integer :: IJK, nc
135 real(dp) :: rr(NDIM)
136 real(dp) :: ellips_fac(NDIM)
137
138 nc = box%n_cell
139
140 ! Create ellipsoidal shape
141 ellips_fac(2:) = 3.0_dp
142 ellips_fac(1) = 1.0_dp
143
144 do kji_do(0,nc+1)
145 rr = af_r_cc(box, [ijk])
146
147 ! Change epsilon in part of the domain
148 if (norm2((rr - 0.5_dp) * ellips_fac) < 0.25_dp) then
149 box%cc(ijk, i_eps) = epsilon_high
150 else
151 box%cc(ijk, i_eps) = 1.0_dp
152 end if
153
154 box%cc(ijk, i_rhs) = 0.0d0
155 box%cc(ijk, i_phi) = 0.0d0
156 end do; close_do
157
158 end subroutine set_init_cond
159
160 ! This routine sets boundary conditions for a box
161 subroutine sides_bc(box, nb, iv, coords, bc_val, bc_type)
162 type(box_t), intent(in) :: box
163 integer, intent(in) :: nb
164 integer, intent(in) :: iv
165 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
166 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
167 integer, intent(out) :: bc_type
168 integer :: nc
169
170 nc = box%n_cell
171
172 ! Below the solution is specified in the approriate ghost cells
173 select case (nb)
174 case (af_neighb_lowx) ! Lower-x direction
175 bc_type = af_bc_dirichlet
176 bc_val = 0.0_dp
177 case (af_neighb_highx) ! Higher-x direction
178 bc_type = af_bc_dirichlet
179 bc_val = 1.0_dp
180 case default
181 bc_type = af_bc_neumann
182 bc_val = 0.0_dp
183 end select
184 end subroutine sides_bc
185
186end program dielectric_test
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