Afivo 0.3
All Classes Namespaces Functions Variables Pages
compressible_flow_wall.f90

Example with compressible flow interacting with a solid wall.

Example with compressible flow interacting with a solid wall

1#include "../src/cpp_macros.h"
2!> \example compressible_flow_wall.f90
3!>
4!> Example with compressible flow interacting with a solid wall
5program compressible_flow_wall
6 use m_af_all
7 implicit none
8
9 integer, parameter :: n_vars = 2+ndim
10 integer, parameter :: box_size = 8
11
12 ! Indices defining the order of the flux variables
13 integer, parameter :: i_rho = 1
14 integer, parameter :: i_mom(NDIM) = [2, 3]
15 integer, parameter :: i_e = 4
16 integer :: i_lsf
17
18 real(dp), parameter :: k_boltzmann = 1.380649e-23_dp ! J/K
19 real(dp), parameter :: dalton = 1.66053906660e-27_dp ! kg
20
21 real(dp) :: euler_gamma = 1.4_dp
22 real(dp) :: cfl_number = 0.5_dp
23 integer :: time_integrator = af_heuns_method
24 real(dp) :: rho_0, u_0, p_0, T_0, sound_speed, mach_number
25 real(dp) :: molecular_mass, domain_size(NDIM)
26 real(dp) :: wedge_angle, wedge_normal(2)
27 integer :: variables(n_vars)
28 integer :: cc_rho, cc_mom(2), cc_e
29 integer :: fluxes(n_vars)
30 integer :: grid_size(NDIM)
31 integer :: iteration
32 real(dp) :: u_init(1, n_vars)
33 type(af_t) :: tree
34 real(dp) :: dt, dt_lim, time, end_time, dt_output
35 character(len=100) :: fname
36 integer :: n, output_cnt
37 logical :: use_amr
38
39 character(len=10) :: var_names(n_vars)
40
41 ! This test case with supersonic flow past a wedge is described in
42 ! https://doi.org/10.3929/ethz-a-005575219 (master thesis)
43
44 domain_size = [0.45_dp, 0.3_dp] ! m
45 grid_size = [12, 8] * box_size
46 end_time = 3e-3_dp ! s
47 euler_gamma = 1.4_dp
48 molecular_mass = 28.9_dp ! daltons
49 p_0 = 101.35e3_dp ! Pa
50 t_0 = 288.9_dp ! K
51 wedge_angle = 15.0_dp * acos(-1.0_dp) / 180.0_dp
52 wedge_normal = [-sin(wedge_angle), cos(wedge_angle)]
53
54 ! gas density in kg/m3
55 rho_0 = p_0 * molecular_mass * dalton / (k_boltzmann * t_0)
56 sound_speed = sqrt(euler_gamma * p_0 / rho_0) ! m/s
57 mach_number = 2.5_dp
58 u_0 = mach_number * sound_speed ! m/s
59
60 u_init(1, :) = [rho_0, u_0, 0.0_dp, p_0]
61 call to_conservative(1, n_vars, u_init)
62
63 var_names(i_rho) = "rho"
64 var_names(i_mom) = ["mom_x", "mom_y"]
65 var_names(i_e) = "E"
66
67 do n = 1, n_vars
68 call af_add_cc_variable(tree, var_names(n), ix=variables(n), n_copies=2)
69 call af_add_fc_variable(tree, "flux", ix=fluxes(n))
70 end do
71
72 call af_set_cc_methods(tree, variables(i_rho), bc_gas_rho)
73 call af_set_cc_methods(tree, variables(i_mom(1)), bc_gas_mom_x)
74 call af_set_cc_methods(tree, variables(i_mom(2)), bc_gas_mom_y)
75 call af_set_cc_methods(tree, variables(i_e), bc_gas_e)
76
77 call af_add_cc_variable(tree, "lsf", ix=i_lsf)
78 call af_set_cc_methods(tree, i_lsf, funcval=set_lsf_box)
79
80 cc_rho = variables(i_rho)
81 cc_mom = variables(i_mom)
82 cc_e = variables(i_e)
83
84 use_amr = .false.
85
86 call af_init(tree, box_size, domain_size, grid_size)
87 call af_loop_box(tree, set_init_conds)
88 call af_gc_tree(tree, variables)
89
90 time = 0.0_dp
91 dt = 0.0_dp ! Start from zero time step
92 output_cnt = 0
93 iteration = 0
94 dt_output = end_time / 100
95
96 do
97 iteration = iteration + 1
98 if (output_cnt * dt_output <= time) then
99 write(fname, "(A,I0)") "output/compressible_flow_wall_" &
100 // dimname // "_", output_cnt
101 call af_write_silo(tree, trim(fname), output_cnt, time, &
102 add_vars = write_primitives, add_names=["xVel","yVel","pres"])
103 output_cnt = output_cnt + 1
104 end if
105
106 call af_advance(tree, dt, dt_lim, time, variables, &
107 time_integrator, forward_euler)
108 dt = cfl_number * dt_lim
109 if (time > end_time) exit
110 end do
111
112 print *, "Number of time steps: ", iteration
113 print *, "Number of grid cells: ", af_num_leaves_used(tree) * box_size**ndim
114 call af_destroy(tree)
115
116contains
117
118 subroutine forward_euler(tree, dt, dt_stiff, dt_lim, time, s_deriv, n_prev, &
119 s_prev, w_prev, s_out, i_step, n_steps)
120 type(af_t), intent(inout) :: tree
121 real(dp), intent(in) :: dt
122 real(dp), intent(in) :: dt_stiff !< Time step for stiff terms
123 real(dp), intent(inout) :: dt_lim
124 real(dp), intent(in) :: time
125 integer, intent(in) :: s_deriv
126 integer, intent(in) :: n_prev !< Number of previous states
127 integer, intent(in) :: s_prev(n_prev) !< Previous states
128 real(dp), intent(in) :: w_prev(n_prev) !< Weights of previous states
129 integer, intent(in) :: s_out
130 integer, intent(in) :: i_step, n_steps
131 real(dp) :: dummy_dt(0)
132
133 call flux_generic_tree(tree, n_vars, variables, s_deriv, fluxes, dt_lim, &
134 max_wavespeed, get_fluxes, flux_dummy_modify, line_modify, &
135 to_primitive, to_conservative, af_limiter_vanleer_t)
136 call flux_update_densities(tree, dt, n_vars, variables, n_vars, variables, fluxes, &
137 s_deriv, n_prev, s_prev, w_prev, s_out, flux_dummy_source, 0, dummy_dt, &
138 set_box_mask)
139 end subroutine forward_euler
140
141 !> Set a mask to true where the solution should be updated
142 subroutine set_box_mask(box, mask)
143 type(box_t), intent(in) :: box
144 logical, intent(out) :: mask(DTIMES(box%n_cell))
145 integer :: nc
146
147 nc = box%n_cell
148 mask = (box%cc(dtimes(1:nc), i_lsf) > 0.0_dp)
149 end subroutine set_box_mask
150
151 subroutine set_init_conds(box)
152 type(box_t), intent(inout) :: box
153 integer :: IJK, nc
154
155 nc = box%n_cell
156 do kji_do(0, nc+1)
157 box%cc(ijk, variables) = u_init(1, :)
158 if (box%cc(ijk, i_lsf) <= 0) then
159 box%cc(ijk, i_mom) = 0.0_dp
160 end if
161 end do; close_do
162 end subroutine set_init_conds
163
164 !> Convert variables to primitive form
165 subroutine to_primitive(n_values, n_vars, u)
166 integer, intent(in) :: n_values, n_vars
167 real(dp), intent(inout) :: u(n_values, n_vars)
168
169 u(:, i_mom(1)) = u(:, i_mom(1))/u(:, i_rho)
170 u(:, i_mom(2)) = u(:, i_mom(2))/u(:, i_rho)
171 u(:, i_e) = (euler_gamma-1.0_dp) * (u(:, i_e) - &
172 0.5_dp*u(:, i_rho)* sum(u(:, i_mom(:))**2, dim=2))
173 end subroutine to_primitive
174
175 !> Convert variables to conservative form
176 subroutine to_conservative(n_values, n_vars, u)
177 integer, intent(in) :: n_values, n_vars
178 real(dp), intent(inout) :: u(n_values, n_vars)
179 real(dp) :: kin_en(n_values)
180 real(dp) :: inv_fac
181 integer :: i
182
183 ! Compute kinetic energy (0.5 * rho * velocity^2)
184 kin_en = 0.5_dp * u(:, i_rho) * sum(u(:, i_mom(:))**2, dim=2)
185
186 ! Compute energy from pressure and kinetic energy
187 inv_fac = 1/(euler_gamma - 1.0_dp)
188 u(:, i_e) = u(:, i_e) * inv_fac + kin_en
189
190 ! Compute momentum from density and velocity components
191 do i = 1, ndim
192 u(:, i_mom(i)) = u(:, i_rho) * u(:, i_mom(i))
193 end do
194 end subroutine to_conservative
195
196 !> Get maximal wave speed
197 subroutine max_wavespeed(n_values, n_var, flux_dim, u, w)
198 integer, intent(in) :: n_values !< Number of cell faces
199 integer, intent(in) :: n_var !< Number of variables
200 integer, intent(in) :: flux_dim !< In which dimension fluxes are computed
201 real(dp), intent(in) :: u(n_values, n_var) !< Primitive variables
202 real(dp), intent(out) :: w(n_values) !< Maximum speed
203 real(dp) :: sound_speeds(n_values)
204
205 sound_speeds = sqrt(euler_gamma * u(:, i_e) / u(:, i_rho))
206 w = sound_speeds + abs(u(:, i_mom(flux_dim)))
207 end subroutine max_wavespeed
208
209 !> Compute fluxes from variables in primitive form
210 subroutine get_fluxes(n_values, n_var, flux_dim, u, flux, box, line_ix, s_deriv)
211 integer, intent(in) :: n_values !< Number of cell faces
212 integer, intent(in) :: n_var !< Number of variables
213 integer, intent(in) :: flux_dim !< In which dimension fluxes are computed
214 real(dp), intent(in) :: u(n_values, n_var)
215 real(dp), intent(out) :: flux(n_values, n_var)
216 type(box_t), intent(in) :: box
217 integer, intent(in) :: line_ix(NDIM-1)
218 integer, intent(in) :: s_deriv
219 real(dp) :: E(n_values), inv_fac
220 integer :: i
221
222 ! Compute left and right flux for conservative variables from the primitive
223 ! reconstructed values.
224
225 ! Density flux
226 flux(:, i_rho) = u(:, i_rho) * u(:, i_mom(flux_dim))
227
228 ! Momentum flux
229 do i = 1, ndim
230 flux(:, i_mom(i)) = u(:, i_rho) * &
231 u(:, i_mom(i)) * u(:, i_mom(flux_dim))
232 end do
233
234 ! Add pressure term
235 flux(:, i_mom(flux_dim)) = flux(:, i_mom(flux_dim)) + u(:, i_e)
236
237 ! Compute energy
238 inv_fac = 1/(euler_gamma-1.0_dp)
239 e = u(:, i_e) * inv_fac + 0.5_dp * u(:, i_rho) * &
240 sum(u(:, i_mom(:))**2, dim=2)
241
242 ! Energy flux
243 flux(:, i_e) = u(:, i_mom(flux_dim)) * (e + u(:, i_e))
244
245 end subroutine get_fluxes
246
247 !> Modify line data to approximate a slip boundary condition
248 subroutine line_modify(n_cc, n_var, cc_line, flux_dim, box, line_ix, s_deriv)
249 integer, intent(in) :: n_cc !< Number of cell centers
250 integer, intent(in) :: n_var !< Number of variables
251 real(dp), intent(inout) :: cc_line(-1:n_cc-2, n_var) !< Line values to modify
252 integer, intent(in) :: flux_dim !< In which dimension fluxes are computed
253 type(box_t), intent(in) :: box !< Current box
254 integer, intent(in) :: line_ix(NDIM-1) !< Index of line for dim /= flux_dim
255 integer, intent(in) :: s_deriv !< State to compute derivatives from
256
257 real(dp) :: lsf(0:box%n_cell+1)
258 real(dp) :: mom(2), mom_par(2), mom_perp(2)
259 integer :: i
260
261 ! Get level set function along the line of the flux computation
262 call flux_get_line_1cc(box, i_lsf, flux_dim, line_ix, lsf)
263
264 if (all(lsf > 0)) return ! no boundary
265
266 do i = 0, box%n_cell
267 if (lsf(i) * lsf(i+1) <= 0) then
268 ! There is an interface
269 if (lsf(i) > 0) then
270 cc_line(i+1, :) = cc_line(i, :)
271 cc_line(i+2, :) = cc_line(i-1, :)
272
273 ! Project momentum onto a component parallel and perpendicular to
274 ! the boundary. The parallel momentum is copied, but for the
275 ! perpendicular moment a negative sign is used so that the
276 ! "average" perpendicular momentum is zero.
277 ! TODO: use gradient of lsf, take distance into account?
278 mom = cc_line(i, i_mom)
279 mom_perp = dot_product(mom, wedge_normal) * wedge_normal
280 mom_par = mom - mom_perp
281 cc_line(i+1, i_mom) = mom_par - mom_perp
282
283 mom = cc_line(i-1, i_mom)
284 mom_perp = dot_product(mom, wedge_normal) * wedge_normal
285 mom_par = mom - mom_perp
286 cc_line(i+2, i_mom) = mom_par - mom_perp
287 else
288 cc_line(i, :) = cc_line(i+1, :)
289 cc_line(i-1, :) = cc_line(i+2, :)
290
291 mom = cc_line(i+1, i_mom)
292 mom_perp = dot_product(mom, wedge_normal) * wedge_normal
293 mom_par = mom - mom_perp
294 cc_line(i, i_mom) = mom_par - mom_perp
295
296 mom = cc_line(i+2, i_mom)
297 mom_perp = dot_product(mom, wedge_normal) * wedge_normal
298 mom_par = mom - mom_perp
299 cc_line(i-1, i_mom) = mom_par - mom_perp
300 end if
301 end if
302 end do
303
304 end subroutine line_modify
305
306 !> Write primitive variables to output
307 subroutine write_primitives(box, new_vars, n_var)
308 type(box_t), intent(in) :: box
309 integer, intent(in) :: n_var
310 real(dp) :: new_vars(DTIMES(0:box%n_cell+1),n_var)
311
312 ! X Velocity
313 new_vars(dtimes(:), 1) = box%cc(dtimes(:), cc_mom(1))/box%cc(dtimes(:), cc_rho)
314 ! Y Velocity
315 new_vars(dtimes(:), 2) = box%cc(dtimes(:), cc_mom(2))/box%cc(dtimes(:), cc_rho)
316 ! Pressure
317 new_vars(dtimes(:), 3) = (euler_gamma-1.0_dp)*(box%cc(dtimes(:), cc_e) - &
318 sum(box%cc(dtimes(:), cc_mom(:))**2, dim=ndim+1) &
319 /(2.0_dp*box%cc(dtimes(:), cc_rho)))
320 end subroutine write_primitives
321
322 !> Boundary conditions for gas density
323 subroutine bc_gas_rho(box, nb, iv, coords, bc_val, bc_type)
324 type(box_t), intent(in) :: box
325 integer, intent(in) :: nb
326 integer, intent(in) :: iv
327 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
328 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
329 integer, intent(out) :: bc_type
330
331 if (nb == af_neighb_lowx) then
332 bc_type = af_bc_dirichlet
333 bc_val = u_init(1, i_rho)
334 else
335 call af_bc_neumann_zero(box, nb, iv, coords, bc_val, bc_type)
336 end if
337 end subroutine bc_gas_rho
338
339 !> Boundary conditions for x-momentum
340 subroutine bc_gas_mom_x(box, nb, iv, coords, bc_val, bc_type)
341 type(box_t), intent(in) :: box
342 integer, intent(in) :: nb
343 integer, intent(in) :: iv
344 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
345 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
346 integer, intent(out) :: bc_type
347
348 if (nb == af_neighb_lowx) then
349 bc_type = af_bc_dirichlet
350 bc_val = u_init(1, i_mom(1))
351 else
352 call af_bc_neumann_zero(box, nb, iv, coords, bc_val, bc_type)
353 end if
354 end subroutine bc_gas_mom_x
355
356 !> Boundary conditions for y-momentum
357 subroutine bc_gas_mom_y(box, nb, iv, coords, bc_val, bc_type)
358 type(box_t), intent(in) :: box
359 integer, intent(in) :: nb
360 integer, intent(in) :: iv
361 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
362 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
363 integer, intent(out) :: bc_type
364
365 if (nb == af_neighb_lowx) then
366 bc_type = af_bc_dirichlet
367 bc_val = u_init(1, i_mom(2))
368 else if (nb == af_neighb_lowy) then
369 bc_type = af_bc_dirichlet
370 bc_val = 0
371 else if (nb == af_neighb_highy) then
372 bc_type = af_bc_dirichlet
373 bc_val = 0
374 else
375 call af_bc_neumann_zero(box, nb, iv, coords, bc_val, bc_type)
376 end if
377 end subroutine bc_gas_mom_y
378
379 !> Boundary conditions for energy
380 subroutine bc_gas_e(box, nb, iv, coords, bc_val, bc_type)
381 type(box_t), intent(in) :: box
382 integer, intent(in) :: nb
383 integer, intent(in) :: iv
384 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
385 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
386 integer, intent(out) :: bc_type
387
388 if (nb == af_neighb_lowx) then
389 bc_type = af_bc_dirichlet
390 bc_val = u_init(1, i_e)
391 else
392 call af_bc_neumann_zero(box, nb, iv, coords, bc_val, bc_type)
393 end if
394 end subroutine bc_gas_e
395
396 subroutine set_lsf_box(box, iv)
397 type(box_t), intent(inout) :: box
398 integer, intent(in) :: iv
399 integer :: IJK, nc
400 real(dp) :: rr(NDIM), norm_dr
401
402 nc = box%n_cell
403 norm_dr = norm2(box%dr)
404
405 do kji_do(0,nc+1)
406 rr = af_r_cc(box, [ijk])
407 box%cc(ijk, iv) = get_lsf(rr)
408 end do; close_do
409 end subroutine set_lsf_box
410
411 !> Level set function
412 real(dp) function get_lsf(rr)
413 real(dp), intent(in) :: rr(NDIM)
414
415 ! Distance from a point to a line
416 ! https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line
417 get_lsf = cos(wedge_angle) * (rr(2) - 0.0_dp) - &
418 sin(wedge_angle) * (rr(1) - 0.15_dp)
419
420 end function get_lsf
421
422end program
Module which contains all Afivo modules, so that a user does not have to include them separately.
Definition: m_af_all.f90:3