An example to test the accuracy of time integration, solving a simple equation: d/dt y = -y.
An example to test the accuracy of time integration, solving a simple equation: d/dt y = -y
1#include "../src/cpp_macros.h"
2
3
4
5
6program time_integration
8
9 implicit none
10
11 integer :: i_phi, i_err
12 integer :: n, n_steps, integrator
13 real(dp), parameter :: initial_value = exp(1.0_dp)
14 real(dp) :: dt, time, end_time
15 real(dp) :: max_error
16 character(len=100) :: argument
17 type(af_t) :: tree
18
19 end_time = 1.0_dp
20 dt = 0.1_dp
21
22
23 if (command_argument_count() > 0) then
24 call get_command_argument(1, argument)
25 read(argument, *) dt
26 end if
27
28 n_steps = nint(end_time/dt)
29 dt = end_time / n_steps
30
31 write(*, "(A20,2A12)") "integrator ", "dt", "max_error"
32
33 do integrator = 1, af_num_integrators
34 time = 0.0_dp
35 call integrate(tree, integrator, time, dt, n_steps)
36 call af_loop_box_arg(tree, set_error, [time])
37 call af_tree_maxabs_cc(tree, i_err, max_error)
38 write(*, "(A20,2E12.4)") af_integrator_names(integrator), dt, max_error
39 call af_destroy(tree)
40 end do
41
42contains
43
44 subroutine set_initial_condition(box)
45 type(box_t), intent(inout) :: box
46 box%cc(dtimes(:), i_phi) = initial_value
47 end subroutine set_initial_condition
48
49 elemental real(dp) function solution(t)
50 real(dp), intent(in) :: t
51 solution = initial_value * exp(-t)
52 end function solution
53
54 subroutine set_error(box, time)
55 type(box_t), intent(inout) :: box
56 real(dp), intent(in) :: time(:)
57 integer :: IJK, nc
58
59 nc = box%n_cell
60 do kji_do(1,nc)
61 box%cc(ijk, i_err) = box%cc(ijk, i_phi) - solution(time(1))
62 end do; close_do
63 end subroutine set_error
64
65 subroutine forward_euler(tree, dt, dt_stiff, dt_lim, time, s_deriv, n_prev, s_prev, &
66 w_prev, s_out, i_step, n_steps)
67 type(af_t), intent(inout) :: tree
68 real(dp), intent(in) :: dt
69 real(dp), intent(in) :: dt_stiff
70 real(dp), intent(inout) :: dt_lim
71 real(dp), intent(in) :: time
72 integer, intent(in) :: s_deriv
73 integer, intent(in) :: n_prev
74 integer, intent(in) :: s_prev(n_prev)
75 real(dp), intent(in) :: w_prev(n_prev)
76 integer, intent(in) :: s_out
77 integer, intent(in) :: i_step
78 integer, intent(in) :: n_steps
79 integer :: lvl, IJK, n, id, nc
80
81 nc = tree%n_cell
82
83
84 do lvl = 1, tree%highest_lvl
85
86 do n = 1, size(tree%lvls(lvl)%leaves)
87 id = tree%lvls(lvl)%leaves(n)
88 associate(cc => tree%boxes(id)%cc)
89 do kji_do(1,nc)
90 cc(ijk, i_phi+s_out) = sum(w_prev * cc(ijk, i_phi+s_prev)) - &
91 dt_stiff * cc(ijk, i_phi+s_deriv)
92 end do; close_do
93 end associate
94 end do
95
96 end do
97
98
99 call af_gc_tree(tree, [i_phi])
100 end subroutine forward_euler
101
102 subroutine integrate(tree, integrator, time, dt, n_steps)
103 type(af_t), intent(inout) :: tree
104 integer, intent(in) :: integrator
105 real(dp), intent(inout) :: time
106 real(dp), intent(in) :: dt
107 integer, intent(in) :: n_steps
108 real(dp) :: dt_lim
109 integer :: n_copies
110 integer, parameter :: box_size = 8
111
112 n_copies = af_advance_num_steps(integrator)
113 call af_add_cc_variable(tree, "phi", ix=i_phi, n_copies=n_copies)
114 call af_add_cc_variable(tree, "err", ix=i_err)
115 call af_set_cc_methods(tree, i_phi, af_bc_neumann_zero)
116 call af_init(tree, box_size, [dtimes(1.0_dp)], [dtimes(box_size)], &
117 periodic=[dtimes(.true.)], mem_limit_gb=1e-3_dp)
118
119 call af_loop_box(tree, set_initial_condition)
120
121 do n = 1, n_steps
122 call af_advance(tree, dt, dt_lim, time, [i_phi], integrator, &
123 forward_euler, implicit_solver)
124 end do
125 end subroutine integrate
126
127
128
129 subroutine implicit_solver(tree, dt_stiff, time, n_prev, s_prev, w_prev, s_out)
130 type(af_t), intent(inout) :: tree
131 real(dp), intent(in) :: dt_stiff
132 real(dp), intent(in) :: time
133 integer, intent(in) :: n_prev
134 integer, intent(in) :: s_prev(n_prev)
135 real(dp), intent(in) :: w_prev(n_prev)
136 integer, intent(in) :: s_out
137
138 integer :: lvl, IJK, n, id, nc
139
140 nc = tree%n_cell
141
142
143 do lvl = 1, tree%highest_lvl
144
145 do n = 1, size(tree%lvls(lvl)%leaves)
146 id = tree%lvls(lvl)%leaves(n)
147 associate(cc => tree%boxes(id)%cc)
148 do kji_do(1,nc)
149 cc(ijk, i_phi+s_out) = sum(w_prev * cc(ijk, i_phi+s_prev)) / &
150 (1 + dt_stiff)
151 end do; close_do
152 end associate
153 end do
154
155 end do
156
157
158 call af_gc_tree(tree, [i_phi])
159 end subroutine implicit_solver
160
161end program time_integration
Module which contains all Afivo modules, so that a user does not have to include them separately.