Example showing how to use different types of boundary conditions.
1#include "../src/cpp_macros.h"
2
3
4
5program boundary_conditions
7
8 implicit none
9
10 integer, parameter :: box_size = 8
11 integer, parameter :: n_iterations = 20
12 integer :: i_phi, i_rho
13
14 type(af_t) :: tree
15 integer :: iter
16 character(len=100) :: fname
17
18 print *, "Running boundary_conditions_" // dimname // ""
19 print *, "Number of threads", af_get_max_threads()
20
21 call af_add_cc_variable(tree, "phi", ix=i_phi)
22 call af_add_cc_variable(tree, "rho", ix=i_rho)
23 call af_set_cc_methods(tree, i_phi, boundary_method)
24 call af_set_cc_methods(tree, i_rho, bc_custom=custom_boundary_method)
25
26
27 call af_init(tree, &
28 box_size, &
29 [dtimes(1.0_dp)], &
30 [dtimes(box_size)])
31
32 call af_print_info(tree)
33
34 do iter = 1, n_iterations
35 if (iter == 1) then
36
37 call af_loop_box(tree, set_initial_cond)
38 else
39
40 call af_loop_box(tree, average_vars)
41 end if
42
43 call af_gc_tree(tree, [i_phi, i_rho])
44
45 write(fname, "(A,I0)") "output/boundary_conditions_" // dimname // "_", iter
46 call af_write_silo(tree, trim(fname))
47 end do
48
49contains
50
51
52 subroutine set_initial_cond(box)
53 type(box_t), intent(inout) :: box
54
55 box%cc(dtimes(:), i_phi) = 0.0_dp
56 box%cc(dtimes(:), i_rho) = 1.0_dp
57 end subroutine set_initial_cond
58
59 subroutine average_vars(box)
60 type(box_t), intent(inout) :: box
61 integer :: IJK, nc, ivs(2)
62 real(dp) :: tmp(DTIMES(box%n_cell), 2)
63
64 nc = box%n_cell
65 ivs = [i_phi, i_rho]
66
67 do kji_do(1,nc)
68#if NDIM == 1
69 tmp(i, :) = 0.5_dp * ( &
70 box%cc(i+1, ivs) + &
71 box%cc(i-1, ivs))
72#elif NDIM == 2
73 tmp(i, j, :) = 0.25_dp * ( &
74 box%cc(i+1, j, ivs) + &
75 box%cc(i-1, j, ivs) + &
76 box%cc(i, j+1, ivs) + &
77 box%cc(i, j-1, ivs))
78#elif NDIM == 3
79 tmp(i, j, k, :) = 1/6.0_dp * ( &
80 box%cc(i+1, j, k, ivs) + &
81 box%cc(i-1, j, k, ivs) + &
82 box%cc(i, j+1, k, ivs) + &
83 box%cc(i, j-1, k, ivs) + &
84 box%cc(i, j, k+1, ivs) + &
85 box%cc(i, j, k-1, ivs))
86#endif
87 end do; close_do
88
89
90 box%cc(dtimes(1:nc), ivs) = 0.5_dp * (&
91 tmp(dtimes(:), :) + box%cc(dtimes(1:nc), ivs))
92 end subroutine average_vars
93
94
95 subroutine boundary_method(box, nb, iv, coords, bc_val, bc_type)
96 type(box_t), intent(in) :: box
97 integer, intent(in) :: nb
98 integer, intent(in) :: iv
99 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
100 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
101 integer, intent(out) :: bc_type
102 integer :: nc
103
104 nc = box%n_cell
105
106
107 select case (nb)
108 case (af_neighb_lowx)
109 bc_type = af_bc_dirichlet
110 bc_val = 1.0_dp
111 case default
112 bc_type = af_bc_neumann
113 bc_val = 0.0_dp
114 end select
115
116 end subroutine boundary_method
117
118
119
120 subroutine custom_boundary_method(box, nb, iv, n_gc, cc)
121 type(box_t), intent(inout) :: box
122 integer, intent(in) :: nb
123 integer, intent(in) :: iv
124 integer, intent(in) :: n_gc
125
126 real(dp), intent(inout), optional :: &
127 cc(DTIMES(1-n_gc:box%n_cell+n_gc))
128
129 integer :: lo(NDIM), hi(NDIM)
130
131 if (n_gc /= 1) error stop "not implemented"
132
133
134 call af_get_index_bc_outside(nb, box%n_cell, 1, lo, hi)
135
136
137 box%cc(dslice(lo,hi), iv) = 0.0_dp
138 end subroutine custom_boundary_method
139
140end program
Module which contains all Afivo modules, so that a user does not have to include them separately.