This example shows how to define and apply numerical stencils.
1#include "../src/cpp_macros.h"
2
3
4
5program stencil_test
7
8 implicit none
9
10 integer, parameter :: box_size = 8
11 real(dp), parameter :: domain_len = 2 * acos(-1.0_dp)
12 integer :: i_phi
13 integer :: i_phi_old
14 type(af_t) :: tree
15 type(ref_info_t) :: refine_info
16 integer :: it
17 character(len=100) :: fname
18
19 integer :: stencil_key = 1
20
21 call af_add_cc_variable(tree, "phi", ix=i_phi)
22 call af_add_cc_variable(tree, "phi_old", ix=i_phi_old)
23 call af_set_cc_methods(tree, i_phi, af_bc_neumann_zero)
24 call af_set_cc_methods(tree, i_phi_old, af_bc_neumann_zero)
25
26
27 call af_init(tree, &
28 box_size, &
29 [dtimes(domain_len)], &
30 [dtimes(box_size)], &
31 periodic=[dtimes(.true.)])
32
33 do
34 call af_loop_box(tree, set_initial_condition)
35 call af_adjust_refinement(tree, refine_routine, refine_info)
36 if (refine_info%n_add == 0) exit
37 end do
38
39 call af_stencil_store(tree, stencil_key, set_stencil)
40
41 do it = 1, 10
42 write(fname, "(A,I0)") "output/stencil_test_" // dimname // "_", it
43
44
45
46 call af_write_silo(tree, trim(fname), it)
47 call af_tree_copy_cc(tree, i_phi, i_phi_old)
48 call af_stencil_apply(tree, stencil_key, i_phi_old, i_phi)
49 call af_gc_tree(tree, [i_phi])
50 end do
51
52contains
53
54
55 subroutine refine_routine(box, cell_flags)
56 type(box_t), intent(in) :: box
57 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
58
59 if (maxval(box%dr) > 2e-2_dp * domain_len) then
60 cell_flags = af_do_ref
61 else
62 cell_flags = af_keep_ref
63 end if
64 end subroutine refine_routine
65
66
67 subroutine set_initial_condition(box)
68 type(box_t), intent(inout) :: box
69 integer :: IJK, nc
70 real(dp) :: rr(NDIM)
71
72 nc = box%n_cell
73
74 do kji_do(0,nc+1)
75 rr = af_r_cc(box, [ijk])
76 box%cc(ijk, i_phi) = product(sin(rr))
77 end do; close_do
78 end subroutine set_initial_condition
79
80 subroutine set_stencil(box, stencil)
81 type(box_t), intent(in) :: box
82 type(stencil_t), intent(inout) :: stencil
83
84 stencil%shape = af_stencil_357
85 stencil%stype = stencil_constant
86 allocate(stencil%c(2*ndim+1))
87
88 stencil%c(1) = 0.0_dp
89 stencil%c(2:) = 1.0_dp / (2 * ndim)
90 end subroutine set_stencil
91
92end program stencil_test
Module which contains all Afivo modules, so that a user does not have to include them separately.