An implicit diffusion example, showing how the multigrid methods can be used to solve the diffusion equation with a backward Euler scheme.
An implicit diffusion example, showing how the multigrid methods can be used to solve the diffusion equation with a backward Euler scheme.
1#include "../src/cpp_macros.h"
2
3
4
5
6program implicit_diffusion
8
9 implicit none
10
11 integer, parameter :: box_size = 8
12 integer :: i_phi
13 integer :: i_rhs
14 integer :: i_tmp
15 integer :: i_err
16
17 real(dp), parameter :: domain_len = 2 * acos(-1.0_dp)
18 real(dp), parameter :: diffusion_coeff = 1.0_dp
19#if NDIM == 1
20 real(dp), parameter :: solution_modes(NDIM) = [1]
21#elif NDIM == 2
22 real(dp), parameter :: solution_modes(NDIM) = [1, 1]
23#elif NDIM == 3
24 real(dp), parameter :: solution_modes(NDIM) = [1, 1, 0]
25#endif
26
27 type(af_t) :: tree
28 type(mg_t) :: mg
29 type(ref_info_t) :: refine_info
30 integer :: time_steps, output_cnt
31 real(dp) :: time, end_time
32 real(dp), parameter :: dt = 0.1_dp
33 character(len=100) :: fname
34
35 print *, "Running implicit_diffusion_" // dimname // ""
36 print *, "Number of threads", af_get_max_threads()
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, "err", ix=i_err)
42
43
44 call af_init(tree, &
45 box_size, &
46 [dtimes(domain_len)], &
47 [dtimes(box_size)], &
48 periodic=[dtimes(.true.)])
49
50 mg%i_phi = i_phi
51 mg%i_rhs = i_rhs
52 mg%i_tmp = i_tmp
53 mg%sides_bc => af_bc_dirichlet_zero
54 mg%helmholtz_lambda = 1/(diffusion_coeff * dt)
55
56
57
58
59 call mg_init(tree, mg)
60
61 call af_print_info(tree)
62
63
64 do
65
66
67 call af_loop_box(tree, set_initial_condition)
68
69
70
71
72
73
74
75 call af_gc_tree(tree, [i_phi])
76
77
78
79
80
81
82
83 call af_adjust_refinement(tree, &
84 refine_routine, &
85 refine_info)
86
87
88 if (refine_info%n_add == 0) exit
89 end do
90
91 output_cnt = 0
92 time = 0
93 end_time = 2.0_dp
94 time_steps = 0
95 time = 0
96
97
98 do
99 time_steps = time_steps + 1
100
101 output_cnt = output_cnt + 1
102 write(fname, "(A,I0)") "output/implicit_diffusion_" // dimname // "_", output_cnt
103
104
105 call af_loop_box_arg(tree, set_error, [time])
106
107
108
109 call af_write_silo(tree, trim(fname), output_cnt, time)
110
111 if (time > end_time) exit
112
113
114 call af_loop_box(tree, set_rhs)
115
116
117 call mg_fas_fmg(tree, &
118 mg, &
119 set_residual = .true., &
120 have_guess = .true.)
121
122 time = time + dt
123 end do
124
125contains
126
127
128 subroutine refine_routine(box, cell_flags)
129 type(box_t), intent(in) :: box
130 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
131
132 if (maxval(box%dr) > 2e-2_dp * domain_len) then
133 cell_flags = af_do_ref
134 else
135 cell_flags = af_keep_ref
136 end if
137 end subroutine refine_routine
138
139
140 subroutine set_initial_condition(box)
141 type(box_t), intent(inout) :: box
142 integer :: IJK, nc
143 real(dp) :: rr(NDIM)
144
145 nc = box%n_cell
146
147 do kji_do(0,nc+1)
148 rr = af_r_cc(box, [ijk])
149 box%cc(ijk, i_phi) = solution(rr, 0.0_dp)
150 end do; close_do
151 end subroutine set_initial_condition
152
153
154 subroutine set_error(box, time)
155 type(box_t), intent(inout) :: box
156 real(dp), intent(in) :: time(:)
157 integer :: IJK, nc
158 real(dp) :: rr(NDIM)
159
160 nc = box%n_cell
161 do kji_do(1,nc)
162 rr = af_r_cc(box, [ijk])
163 box%cc(ijk, i_err) = &
164 box%cc(ijk, i_phi) - solution(rr, time(1))
165 end do; close_do
166 end subroutine set_error
167
168 function solution(rr, t) result(sol)
169 real(dp), intent(in) :: rr(NDIM), t
170 real(dp) :: sol, tmp(NDIM)
171
172 tmp = solution_modes * rr
173 sol = 1 + product(cos(tmp)) * &
174 exp(-sum(solution_modes**2) * diffusion_coeff * t)
175 end function solution
176
177
178 subroutine set_rhs(box)
179 type(box_t), intent(inout) :: box
180 integer :: nc
181
182 nc = box%n_cell
183 box%cc(dtimes(1:nc), i_rhs) = -1/(dt * diffusion_coeff) * &
184 box%cc(dtimes(1:nc), i_phi)
185 end subroutine set_rhs
186
187end program implicit_diffusion
Module which contains all Afivo modules, so that a user does not have to include them separately.