Example showing how to use multigrid for ensuring a vector field is divergence free.
Example showing how to use multigrid for ensuring a vector field is divergence free
1#include "../src/cpp_macros.h"
2
3
4
5
6program poisson_div_cleaning
9
10 implicit none
11
12 integer, parameter :: box_size = 8
13 integer, parameter :: n_iterations = 10
14 integer :: domain_size(NDIM)
15 real(dp) :: domain_length(NDIM)
16 real(dp) :: domain_origin(NDIM)
17 integer :: i_phi
18 integer :: i_div_old
19 integer :: i_div_new
20 integer :: i_tmp
21 integer :: if_B_old
22 integer :: if_B_new
23 integer :: i_B(3)
24 type(af_t) :: tree
25 type(ref_info_t) :: refine_info
26 integer :: mg_iter
27 character(len=100) :: fname
28 type(mg_t) :: mg
29
30 real(dp), parameter :: alphar_ = 0.64d0
31 real(dp), parameter :: ar_ = 0.25d0
32 real(dp), parameter :: B0 = 1.0_dp
33
34 print *, "Running poisson_div_cleaning_" // dimname
35 print *, "Number of threads", af_get_max_threads()
36
37 call af_add_cc_variable(tree, "phi", ix=i_phi)
38 call af_add_cc_variable(tree, "div_old", ix=i_div_old)
39 call af_add_cc_variable(tree, "tmp", ix=i_tmp)
40 call af_add_cc_variable(tree, "div_new", ix=i_div_new)
41 call af_add_cc_variable(tree, "Bx", ix=i_b(1))
42 call af_add_cc_variable(tree, "By", ix=i_b(2))
43 call af_add_cc_variable(tree, "Bz", ix=i_b(3))
44 call af_add_fc_variable(tree, "B_old", ix=if_b_old)
45 call af_add_fc_variable(tree, "B_new", ix=if_b_new)
46
47 domain_length(:) = 1.0_dp
48 domain_origin(:) = -0.5_dp * domain_length(:)
49 domain_size(:) = box_size
50
51 call af_init(tree, &
52 box_size, &
53 domain_origin+domain_length, &
54 domain_size, &
55 r_min=domain_origin)
56
57 do
58 call af_loop_box(tree, set_initial_condition)
59 call af_adjust_refinement(tree, refine_routine, refine_info)
60 if (refine_info%n_add == 0) exit
61 end do
62
63 call af_print_info(tree)
64
65 mg%i_phi = i_phi
66 mg%i_rhs = i_div_old
67 mg%i_tmp = i_tmp
68 mg%sides_bc => af_bc_dirichlet_zero
69
70 call mg_init(tree, mg)
71
72 do mg_iter = 1, n_iterations
73 call mg_fas_fmg(tree, mg, set_residual=.true., have_guess=(mg_iter>1))
74 call af_loop_box(tree, set_error)
75 write(fname, "(A,I0)") "output/poisson_div_cleaning_" // &
76 dimname // "_", mg_iter
77 call af_write_silo(tree, trim(fname))
78 end do
79
80contains
81
82
83 subroutine refine_routine(box, cell_flags)
84 type(box_t), intent(in) :: box
85 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
86
87 if (box%lvl < 3) then
88 cell_flags = af_do_ref
89 else
90 cell_flags = af_keep_ref
91 end if
92 end subroutine refine_routine
93
94 real(dp) function Bx(xyz)
95 real(dp), intent(in) :: xyz(NDIM)
96 real(dp) :: r, t, phi, Bphi
97 r = norm2(xyz(1:2))
98 t = r / ar_
99 bphi = b0 / (t * (1 + t**2)**alphar_) * &
100 sqrt(((1 + t**2)**(2 * alphar_) - 1 - 2 * alphar_ * t**2) / &
101 (2 * alphar_ - 1))
102 phi = atan2(xyz(2), xyz(1))
103 bx = -sin(phi) * bphi
104 end function bx
105
106 real(dp) function By(xyz)
107 real(dp), intent(in) :: xyz(NDIM)
108 real(dp) :: r, t, phi, Bphi
109 r = norm2(xyz(1:2))
110 t = r / ar_
111 bphi = b0 / (t * (1 + t**2)**alphar_) * &
112 sqrt(((1 + t**2)**(2 * alphar_) - 1 - 2 * alphar_ * t**2) / &
113 (2 * alphar_ - 1))
114 phi = atan2(xyz(2), xyz(1))
115 by = cos(phi) * bphi
116 end function by
117
118 real(dp) function Bz(xyz)
119 real(dp), intent(in) :: xyz(NDIM)
120 real(dp) :: r, t
121 r = norm2(xyz(1:2))
122 t = r / ar_
123 bz = b0 / (1 + t**2)**alphar_
124 end function bz
125
126
127 subroutine set_initial_condition(box)
128 type(box_t), intent(inout) :: box
129 integer :: IJK, nc
130 real(dp) :: xyz(NDIM)
131
132 nc = box%n_cell
133
134
135 do k = 1, nc; do j = 1, nc; do i = 1, nc+1
136 box%fc(ijk, 1, if_b_old) = bx(af_r_fc(box, 1, [ijk]))
137 end do; end do; end do
138
139
140 do k = 1, nc; do j = 1, nc+1; do i = 1, nc
141 box%fc(ijk, 2, if_b_old) = by(af_r_fc(box, 2, [ijk]))
142 end do; end do; end do
143
144
145 do k = 1, nc+1; do j = 1, nc; do i = 1, nc
146 box%fc(ijk, 3, if_b_old) = bz(af_r_fc(box, 3, [ijk]))
147 end do; end do; end do
148
149 do kji_do(1, nc)
150 xyz = af_r_cc(box, [ijk])
151
152
153 box%cc(ijk, i_b(:)) = [bx(xyz), by(xyz), bz(xyz)]
154
155
156 box%cc(ijk, i_div_old) = &
157 (box%fc(i+1, j, k, 1, if_b_old) - &
158 box%fc(i, j, k, 1, if_b_old)) / box%dr(1) + &
159 (box%fc(i, j+1, k, 2, if_b_old) - &
160 box%fc(i, j, k, 2, if_b_old)) / box%dr(2) + &
161 (box%fc(i, j, k+1, 3, if_b_old) - &
162 box%fc(i, j, k, 3, if_b_old)) / box%dr(3)
163 end do; close_do
164 end subroutine set_initial_condition
165
166 subroutine set_error(box)
167 type(box_t), intent(inout) :: box
168 integer :: IJK, nc
169 real(dp) :: inv_dr(NDIM)
170
171 nc = box%n_cell
172 inv_dr = 1 / box%dr
173
174 do k = 1, nc; do j = 1, nc; do i = 1, nc+1
175
176 box%fc(ijk, 1, if_b_new) = box%fc(ijk, 1, if_b_old) - inv_dr(1) * &
177 (box%cc(i, j, k, i_phi) - box%cc(i-1, j, k, i_phi))
178 end do; end do; end do
179
180 do k = 1, nc; do j = 1, nc+1; do i = 1, nc
181
182 box%fc(ijk, 2, if_b_new) = box%fc(ijk, 2, if_b_old) - inv_dr(2) * &
183 (box%cc(i, j, k, i_phi) - box%cc(i, j-1, k, i_phi))
184 end do; end do; end do
185
186 do k = 1, nc+1; do j = 1, nc; do i = 1, nc
187
188 box%fc(ijk, 3, if_b_new) = box%fc(ijk, 3, if_b_old) - inv_dr(3) * &
189 (box%cc(i, j, k, i_phi) - box%cc(i, j, k-1, i_phi))
190 end do; end do; end do
191
192
193 do kji_do(1, nc)
194 box%cc(ijk, i_div_new) = &
195 (box%fc(i+1, j, k, 1, if_b_new) - &
196 box%fc(i, j, k, 1, if_b_new)) / box%dr(1) + &
197 (box%fc(i, j+1, k, 2, if_b_new) - &
198 box%fc(i, j, k, 2, if_b_new)) / box%dr(2) + &
199 (box%fc(i, j, k+1, 3, if_b_new) - &
200 box%fc(i, j, k, 3, if_b_new)) / box%dr(3)
201 end do; close_do
202 end subroutine set_error
203
204end program poisson_div_cleaning
Module which contains all Afivo modules, so that a user does not have to include them separately.
This module can be used to construct solutions consisting of one or more Gaussians.