Afivo 0.3
All Classes Namespaces Functions Variables Pages
implicit_diffusion.f90

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!> \example implicit_diffusion.f90
3!>
4!> An implicit diffusion example, showing how the multigrid methods can be used
5!> to solve the diffusion equation with a backward Euler scheme.
6program implicit_diffusion
7 use m_af_all
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 ! Initialize tree
44 call af_init(tree, & ! Tree to initialize
45 box_size, & ! A box contains box_size**DIM cells
46 [dtimes(domain_len)], &
47 [dtimes(box_size)], &
48 periodic=[dtimes(.true.)])
49
50 mg%i_phi = i_phi ! Solution variable
51 mg%i_rhs = i_rhs ! Right-hand side variable
52 mg%i_tmp = i_tmp ! Variable for temporary space
53 mg%sides_bc => af_bc_dirichlet_zero ! Method for boundary conditions
54 mg%helmholtz_lambda = 1/(diffusion_coeff * dt)
55
56 ! This routine does not initialize the multigrid fields boxes%i_phi,
57 ! boxes%i_rhs and boxes%i_tmp. These fields will be initialized at the
58 ! first call of mg_fas_fmg
59 call mg_init(tree, mg)
60
61 call af_print_info(tree)
62
63 ! Set up the initial conditions
64 do
65
66 ! For each box, set the initial conditions
67 call af_loop_box(tree, set_initial_condition)
68
69 ! Fill ghost cells for variables i_phi on the sides of all boxes, using
70 ! af_gc_interp on refinement boundaries: Interpolation between fine
71 ! points and coarse neighbors to fill ghost cells near refinement
72 ! boundaries.
73 ! Fill ghost cells near physical boundaries using Dirichlet zero
74
75 call af_gc_tree(tree, [i_phi])
76
77 ! Adjust the refinement of a tree using refine_routine (see below) for grid
78 ! refinement.
79 ! Routine af_adjust_refinement sets the bit af_bit_new_children for each box
80 ! that is refined. On input, the tree should be balanced. On output,
81 ! the tree is still balanced, and its refinement is updated (with at most
82 ! one level per call).
83 call af_adjust_refinement(tree, & ! tree
84 refine_routine, & ! Refinement function
85 refine_info) ! Information about refinement
86
87 ! If no new boxes have been added, exit the loop
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 ! Starting simulation
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 ! Call procedure set_error (see below) for each box in tree, with argument time
105 call af_loop_box_arg(tree, set_error, [time])
106
107 ! Write the cell centered data of tree to a file.
108 ! Only the leaves of the tree are used.
109 call af_write_silo(tree, trim(fname), output_cnt, time)
110
111 if (time > end_time) exit
112
113 ! Call set_rhs (see below) for each box in tree
114 call af_loop_box(tree, set_rhs)
115
116 ! Perform a FAS-FMG cycle (full approximation scheme, full multigrid).
117 call mg_fas_fmg(tree, & ! Tree to do multigrid on
118 mg, & ! Multigrid options
119 set_residual = .true., & ! If true, store residual in i_tmp
120 have_guess = .true.) ! If false, start from phi = 0
121
122 time = time + dt
123 end do
124
125contains
126
127 !> Return the refinement flag for box
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 !> This routine sets the initial conditions for each box
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 !> This routine computes the error in i_phi
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 !> This routine computes the right hand side per box
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.
Definition: m_af_all.f90:3