An incompressible flow example.
1#include "../src/cpp_macros.h"
2
3
4
5program incompressible_flow
7
8 implicit none
9
10 integer, parameter :: box_size = 8
11 integer, parameter :: i_mom(NDIM) = [1, 2]
12 integer, parameter :: n_vars = 2
13 integer :: i_pressure = -1
14 real(dp), parameter :: domain_len = 2.0_dp
15 character(len=2), parameter :: mom_names(NDIM) = ["ux", "uy"]
16
17
18 real(dp) :: ux_top = 1.0_dp
19
20
21 real(dp) :: nu = 0.1_dp
22
23
24 real(dp) :: source_ux = 0.0_dp
25
26 type(mg_t) :: mg
27
28 integer, parameter :: integrator = af_midpoint_method
29 integer, parameter :: n_states = af_advance_num_steps(integrator)
30
31 integer :: global_iv_velocity(NDIM)
32 real(dp) :: global_dt
33 integer :: i, it
34 integer :: flow_variables(NDIM)
35 integer :: fluxes(NDIM)
36 type(af_t) :: tree
37 integer :: output_cnt
38 real(dp) :: dt, dt_lim, dt_diff, time
39 real(dp) :: end_time
40 real(dp) :: dt_output
41 character(len=100) :: fname, test_case
42
43 print *, "Running incompressible_flow_" // dimname // ""
44 print *, "Number of threads", af_get_max_threads()
45
46 do i = 1, ndim
47 call af_add_cc_variable(tree, mom_names(i), ix=flow_variables(i), &
48 n_copies=n_states)
49 call af_add_fc_variable(tree, "flux", ix=fluxes(i))
50 end do
51
52 call af_add_cc_variable(tree, "pressure", ix=i_pressure)
53
54 call af_add_cc_variable(tree, "phi", ix=mg%i_phi)
55 call af_add_cc_variable(tree, "rhs", ix=mg%i_rhs)
56 call af_add_cc_variable(tree, "tmp", ix=mg%i_tmp)
57 mg%sides_bc => af_bc_neumann_zero
58 mg%subtract_mean = .false.
59
60 test_case = "channel_flow"
61
62 if (command_argument_count() > 0) then
63 call get_command_argument(1, test_case)
64 end if
65
66 call af_set_cc_methods(tree, i_pressure, af_bc_neumann_zero)
67
68 select case (test_case)
69 case ("cavity_flow")
70 ux_top = 1.0_dp
71 nu = 0.1_dp
72 source_ux = 0.0_dp
73
74 call af_set_cc_methods(tree, flow_variables(1), bc_ux)
75 call af_set_cc_methods(tree, flow_variables(2), af_bc_dirichlet_zero)
76
77 call af_init(tree, box_size, [dtimes(domain_len)], [dtimes(box_size)])
78 case ("channel_flow")
79 nu = 0.1_dp
80 source_ux = 1.0_dp
81
82 call af_set_cc_methods(tree, flow_variables(1), af_bc_dirichlet_zero)
83 call af_set_cc_methods(tree, flow_variables(2), af_bc_dirichlet_zero)
84 call af_init(tree, box_size, [dtimes(domain_len)], [dtimes(box_size)], &
85 periodic=[.true., .false.])
86 case default
87 error stop "Invalid test case, choices: cavity_flow, channel_flow"
88 end select
89
90 call mg_init(tree, mg)
91
92 it = 0
93 output_cnt = 0
94 time = 0
95 dt_output = 1.0_dp
96 end_time = 30.0_dp
97
98 call af_refine_up_to_lvl(tree, 3)
99 call af_loop_box(tree, set_initial_condition)
100 call af_gc_tree(tree, flow_variables)
101
102 dt_diff = 0.5_dp * 0.25_dp * af_min_dr(tree)**2 / nu
103 dt = min(dt_diff, 1e-6_dp)
104 dt_lim = 0.0_dp
105
106
107 do while (time < end_time)
108 if (output_cnt * dt_output <= time) then
109 output_cnt = output_cnt + 1
110 write(fname, "(A,I0)") "output/incompressible_flow_" // dimname // "_", output_cnt
111 call af_gc_tree(tree, flow_variables)
112 call af_write_silo(tree, trim(fname), output_cnt, time)
113 print *, it, dt, dt_lim
114 end if
115
116 call af_advance(tree, dt, dt_lim, time, flow_variables, &
117 af_midpoint_method, forward_euler)
118 dt = 0.9_dp * min(dt_lim, dt_diff, 2.0_dp * dt)
119 it = it + 1
120
121 end do
122
123contains
124
125 subroutine set_initial_condition(box)
126 type(box_t), intent(inout) :: box
127 integer :: IJK, nc
128
129 nc = box%n_cell
130 do kji_do(0,nc+1)
131 box%cc(ijk, flow_variables) = 0.0_dp
132 box%cc(ijk, i_pressure) = 0.0_dp
133 end do; close_do
134 end subroutine set_initial_condition
135
136 subroutine forward_euler(tree, dt, dt_stiff, dt_lim, time, s_deriv, n_prev, &
137 s_prev, w_prev, s_out, i_step, n_steps)
138 type(af_t), intent(inout) :: tree
139 real(dp), intent(in) :: dt
140 real(dp), intent(in) :: dt_stiff
141 real(dp), intent(in) :: time
142 real(dp), intent(inout) :: dt_lim
143 integer, intent(in) :: s_deriv
144 integer, intent(in) :: n_prev
145 integer, intent(in) :: s_prev(n_prev)
146 real(dp), intent(in) :: w_prev(n_prev)
147 integer, intent(in) :: s_out
148 integer, intent(in) :: i_step, n_steps
149 real(dp) :: dummy_dt(0)
150
151 call flux_generic_tree(tree, n_vars, flow_variables, s_deriv, fluxes, dt_lim, &
152 max_wavespeed, get_flux_lr, flux_add_diffusion, flux_dummy_line_modify, &
153 flux_dummy_conversion, flux_dummy_conversion, af_limiter_vanleer_t)
154
155 call flux_update_densities(tree, dt, n_vars, flow_variables, n_vars, &
156 flow_variables, fluxes, s_deriv, n_prev, s_prev, w_prev, s_out, &
157 source_term, 0, dummy_dt)
158
159 call remove_velocity_divergence(tree, flow_variables+s_out, dt)
160
161 end subroutine forward_euler
162
163 subroutine source_term(box, dt, n_vars, i_cc, s_deriv, s_out, n_dt, dt_lim, mask)
164 type(box_t), intent(inout) :: box
165 real(dp), intent(in) :: dt
166 integer, intent(in) :: n_vars
167 integer, intent(in) :: i_cc(n_vars)
168 integer, intent(in) :: s_deriv
169 integer, intent(in) :: s_out
170 integer, intent(in) :: n_dt
171 real(dp), intent(inout) :: dt_lim(n_dt)
172 logical, intent(in) :: mask(DTIMES(box%n_cell))
173
174 integer :: nc
175 nc = box%n_cell
176
177 box%cc(dtimes(1:nc), i_cc(i_mom(1)) + s_out) = &
178 box%cc(dtimes(1:nc), i_cc(i_mom(1)) + s_out) + dt * source_ux
179 end subroutine source_term
180
181 subroutine max_wavespeed(n_values, n_var, flux_dim, u, w)
182 integer, intent(in) :: n_values
183 integer, intent(in) :: n_var
184 integer, intent(in) :: flux_dim
185 real(dp), intent(in) :: u(n_values, n_var)
186 real(dp), intent(out) :: w(n_values)
187
188 w = abs(u(:, i_mom(flux_dim)))
189 end subroutine max_wavespeed
190
191 subroutine get_flux_lr(n_values, n_var, flux_dim, u, flux, box, line_ix, s_deriv)
192 integer, intent(in) :: n_values
193 integer, intent(in) :: n_var
194 integer, intent(in) :: flux_dim
195 real(dp), intent(in) :: u(n_values, n_var)
196 real(dp), intent(out) :: flux(n_values, n_var)
197 type(box_t), intent(in) :: box
198 integer, intent(in) :: line_ix(NDIM-1)
199 integer, intent(in) :: s_deriv
200 integer :: i
201
202 do i = 1, ndim
203
204 flux(:, i_mom(i)) = u(:, i_mom(i)) * u(:, i_mom(flux_dim))
205 end do
206 end subroutine get_flux_lr
207
208 subroutine flux_add_diffusion(n_values, n_var, flux_dim, flux, box, line_ix, s_deriv)
209 integer, intent(in) :: n_values
210 integer, intent(in) :: n_var
211 integer, intent(in) :: flux_dim
212 real(dp), intent(inout) :: flux(n_values, n_var)
213 type(box_t), intent(in) :: box
214 integer, intent(in) :: line_ix(NDIM-1)
215 integer, intent(in) :: s_deriv
216 real(dp) :: inv_dr(NDIM)
217 real(dp) :: cc(0:n_values, 2), u_diff(n_values)
218
219 integer :: i
220
221 inv_dr = 1/box%dr
222
223 call flux_get_line_cc(box, [flow_variables+s_deriv], flux_dim, line_ix, cc)
224
225 do i = 1, ndim
226
227 u_diff = cc(1:, i_mom(i)) - cc(:n_values-1, i_mom(i))
228 flux(:, i_mom(i)) = flux(:, i_mom(i)) - nu * u_diff * inv_dr(i)
229 end do
230 end subroutine flux_add_diffusion
231
232 subroutine remove_velocity_divergence(tree, iv_velocity, dt)
233 type(af_t), intent(inout) :: tree
234 integer, intent(in) :: iv_velocity(NDIM)
235 real(dp), intent(in) :: dt
236 logical :: use_4th_order = .true.
237 real(dp) :: tmp
238 logical, parameter :: print_divergence = .false.
239
240 global_iv_velocity = iv_velocity
241 global_dt = dt
242
243
244 call af_tree_copy_cc(tree, i_pressure, mg%i_phi)
245 call af_loop_box(tree, box_correct_velocity, .true.)
246
247
248 call af_restrict_ref_boundary(tree, iv_velocity)
249
250
251 if (use_4th_order) then
252 call af_loop_tree(tree, box_set_rhs_4th_order, .true.)
253 else
254 call af_loop_tree(tree, box_set_rhs_2nd_order, .true.)
255 end if
256
257 if (print_divergence) then
258 call af_tree_maxabs_cc(tree, mg%i_rhs, tmp)
259 print *, "DIVERGENCE", tmp, iv_velocity(1)
260 end if
261
262
263 call af_tree_clear_cc(tree, mg%i_phi)
264
265
266 call mg_fas_vcycle(tree, mg, .true.)
267
268
269
270 call af_loop_box(tree, box_correct_velocity, .true.)
271
272
273 call af_tree_apply(tree, i_pressure, mg%i_phi, '+')
274
275 end subroutine remove_velocity_divergence
276
277 subroutine box_set_rhs_2nd_order(tree, id)
278 type(af_t), intent(inout) :: tree
279 integer, intent(in) :: id
280 integer :: IJK, nc, iv(NDIM)
281 real(dp) :: tmp(NDIM)
282
283 associate(box => tree%boxes(id))
284 tmp = 0.5_dp / (box%dr * global_dt)
285 iv = global_iv_velocity
286 nc = box%n_cell
287 call af_gc_box(tree, id, iv, .false.)
288
289 do kji_do(1,nc)
290 box%cc(ijk, mg%i_rhs) = &
291 tmp(1) * (box%cc(i+1, j, iv(1)) - box%cc(i-1, j, iv(1))) + &
292 tmp(2) * (box%cc(i, j+1, iv(2)) - box%cc(i, j-1, iv(2)))
293 end do; close_do
294 end associate
295 end subroutine box_set_rhs_2nd_order
296
297 subroutine box_set_rhs_4th_order(tree, id)
298 type(af_t), intent(inout) :: tree
299 integer, intent(in) :: id
300 integer :: IJK, nc, iv(NDIM)
301 real(dp) :: tmp(NDIM)
302 real(dp) :: cc(DTIMES(-1:tree%n_cell+2), n_vars)
303
304 associate(box => tree%boxes(id))
305 nc = box%n_cell
306 tmp = 1/(12 * box%dr * global_dt)
307 iv = global_iv_velocity
308 call af_gc2_box(tree, id, iv, cc)
309
310 do kji_do(1,nc)
311 box%cc(ijk, mg%i_rhs) = &
312 tmp(1) * (8 * (cc(i+1, j, 1) - cc(i-1, j, 1)) - &
313 (cc(i+2, j, 1) - cc(i-2, j, 1))) + &
314 tmp(2) * (8 * (cc(i, j+1, 2) - cc(i, j-1, 2)) - &
315 (cc(i, j+2, 2) - cc(i, j-2, 2)))
316 end do; close_do
317 end associate
318 end subroutine box_set_rhs_4th_order
319
320 subroutine box_correct_velocity(box)
321 type(box_t), intent(inout) :: box
322 integer :: IJK, nc, iv(NDIM)
323 real(dp) :: tmp(NDIM)
324
325 iv = global_iv_velocity
326 tmp = 0.5_dp * global_dt / box%dr
327
328 nc = box%n_cell
329 do kji_do(1,nc)
330 box%cc(ijk, iv(1)) = box%cc(ijk, iv(1)) - tmp(1) * &
331 (box%cc(i+1, j, mg%i_phi) - box%cc(i-1, j, mg%i_phi))
332 box%cc(ijk, iv(2)) = box%cc(ijk, iv(2)) - tmp(2) * &
333 (box%cc(i, j+1, mg%i_phi) - box%cc(i, j-1, mg%i_phi))
334 end do; close_do
335 end subroutine box_correct_velocity
336
337 subroutine bc_ux(box, nb, iv, coords, bc_val, bc_type)
338 type(box_t), intent(in) :: box
339 integer, intent(in) :: nb
340 integer, intent(in) :: iv
341 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
342 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
343 integer, intent(out) :: bc_type
344
345 bc_type = af_bc_dirichlet
346 if (nb == af_neighb_highy) then
347 bc_val = ux_top
348 else
349 bc_val = 0.0_dp
350 end if
351 end subroutine bc_ux
352
353end program incompressible_flow
Module which contains all Afivo modules, so that a user does not have to include them separately.