Example of solving a Helmholtz equation with a stencil that varies in time. Useful ingredient for the implicit solution of diffusion equations.
Example of solving a Helmholtz equation with a stencil that varies in time. Useful ingredient for the implicit solution of diffusion equations.
1#include "../src/cpp_macros.h"
2
3
4
5
6program helmholtz_variable_stencil
8
9 implicit none
10
11 integer, parameter :: box_size = 8
12 integer :: i_err
13
14 real(dp), parameter :: domain_len = 2 * acos(-1.0_dp)
15 real(dp), parameter :: diffusion_coeff = 1.0_dp
16#if NDIM == 1
17 real(dp), parameter :: solution_modes(NDIM) = [1]
18#elif NDIM == 2
19 real(dp), parameter :: solution_modes(NDIM) = [1, 1]
20#elif NDIM == 3
21 real(dp), parameter :: solution_modes(NDIM) = [1, 1, 0]
22#endif
23
24 type(af_t) :: tree
25 type(mg_t) :: mg
26 integer :: i, n, k_factor, n_steps, n_steps_initial
27 real(dp) :: time, dt, dt_initial
28 character(len=100) :: fname
29
30 print *, "Running helmholtz_variable_stencil_" // dimname // ""
31 print *, "Number of threads", af_get_max_threads()
32
33 dt_initial = 0.1_dp
34 n_steps_initial = 10
35
36 call af_add_cc_variable(tree, "phi", ix=mg%i_phi)
37 call af_add_cc_variable(tree, "rhs", ix=mg%i_rhs)
38 call af_add_cc_variable(tree, "tmp", ix=mg%i_tmp)
39 call af_add_cc_variable(tree, "err", ix=i_err)
40
41 call af_init(tree, box_size, [dtimes(domain_len)], &
42 [dtimes(box_size)], periodic=[dtimes(.true.)])
43
44 mg%sides_bc => af_bc_dirichlet_zero
45 mg%helmholtz_lambda = 0.
46 mg%helmholtz_lambda = 1/(diffusion_coeff * dt_initial)
47
48 call mg_init(tree, mg)
49 call af_print_info(tree)
50 call af_refine_up_to_lvl(tree, 3)
51
52 i = 0
53
54 do k_factor = 1, 10
55 time = 0
56 dt = dt_initial / k_factor
57 n_steps = n_steps_initial * k_factor
58 mg%helmholtz_lambda = 1/(diffusion_coeff * dt)
59 print *, k_factor, dt, n_steps
60
61 if (k_factor > 1) then
62 call mg_update_operator_stencil(tree, mg, .false., .true.)
63 end if
64
65 call af_loop_box(tree, set_initial_condition)
66 call af_gc_tree(tree, [mg%i_phi])
67
68 do n = 1, n_steps
69 call af_loop_box(tree, set_rhs)
70 call mg_fas_fmg(tree, mg, .true., .true.)
71 time = time + dt
72 call af_loop_box_arg(tree, set_error, [time])
73 end do
74
75 write(fname, "(A,I0)") "output/helmholtz_variable_stencil_" // dimname // "_", k_factor
76 call af_write_silo(tree, trim(fname), k_factor, time)
77 end do
78
79contains
80
81
82 subroutine set_initial_condition(box)
83 type(box_t), intent(inout) :: box
84 integer :: IJK, nc
85 real(dp) :: rr(NDIM)
86
87 nc = box%n_cell
88
89 do kji_do(0,nc+1)
90 rr = af_r_cc(box, [ijk])
91 box%cc(ijk, mg%i_phi) = solution(rr, 0.0_dp)
92 end do; close_do
93 end subroutine set_initial_condition
94
95
96 subroutine set_error(box, time)
97 type(box_t), intent(inout) :: box
98 real(dp), intent(in) :: time(:)
99 integer :: IJK, nc
100 real(dp) :: rr(NDIM)
101
102 nc = box%n_cell
103 do kji_do(1,nc)
104 rr = af_r_cc(box, [ijk])
105 box%cc(ijk, i_err) = &
106 box%cc(ijk, mg%i_phi) - solution(rr, time(1))
107 end do; close_do
108 end subroutine set_error
109
110 function solution(rr, t) result(sol)
111 real(dp), intent(in) :: rr(NDIM), t
112 real(dp) :: sol, tmp(NDIM)
113
114 tmp = solution_modes * rr
115 sol = 1 + product(cos(tmp)) * &
116 exp(-sum(solution_modes**2) * diffusion_coeff * t)
117 end function solution
118
119
120 subroutine set_rhs(box)
121 type(box_t), intent(inout) :: box
122 integer :: nc
123
124 nc = box%n_cell
125 box%cc(dtimes(1:nc), mg%i_rhs) = -1/(dt * diffusion_coeff) * &
126 box%cc(dtimes(1:nc), mg%i_phi)
127 end subroutine set_rhs
128
129end program helmholtz_variable_stencil
Module which contains all Afivo modules, so that a user does not have to include them separately.
