This example shows how to solve reaction-diffusion equations with different time step methods.
This example shows how to solve reaction-diffusion equations with different time step methods.
1#include "../src/cpp_macros.h"
2
3
4
5
6program reaction_diffusion
8
9 implicit none
10
11 integer, parameter :: box_size = 8
12 integer :: nx_min = 100
13 integer :: i_u
14 integer :: i_v
15 integer :: i_phi1, i_phi2
16 integer :: i_rhs1, i_rhs2
17 integer :: i_tmp
18
19 real(dp) :: domain_len = 1.0_dp
20 type(af_t) :: tree
21 type(mg_t) :: mg1, mg2
22 type(ref_info_t) :: refine_info
23 integer :: it, output_cnt
24 integer :: my_unit
25 real(dp) :: time
26 real(dp) :: end_time = 2.0_dp
27 real(dp) :: dt = 1e-3_dp
28 real(dp) :: dt_output = 1e-2_dp
29 real(dp) :: u_rng_ampl = 0.0_dp
30 real(dp) :: v_rng_ampl = 0.0_dp
31 logical :: periodic(3) = .false.
32 character(len=100) :: fname
33 logical :: file_exists
34
35 character(len=20) :: time_integrator = "imex"
36
37
38
39 character(len=20) :: equation_type = "schnakenberg"
40
41
42
43
44 real(dp) :: alpha = 0.1305_dp
45 real(dp) :: beta = 0.7695_dp
46 real(dp) :: D1 = 0.05_dp
47 real(dp) :: D2 = 1.0_dp
48 real(dp) :: kappa = 100.0_dp
49
50
51 real(dp) :: gs_F = 0.046d0
52 real(dp) :: gs_k = 0.063d0
53
54 namelist /settings/ alpha, beta, d1, d2, kappa, end_time, dt, &
55 dt_output, nx_min, time_integrator, periodic, u_rng_ampl, &
56 v_rng_ampl, equation_type, domain_len, gs_f, gs_k
57
58 print *, "Running reaction_diffusion_" // dimname // ""
59 print *, "Number of threads", af_get_max_threads()
60
61 inquire(file="reaction_diffusion.txt", exist=file_exists)
62 if (file_exists) then
63 open(newunit=my_unit, file="reaction_diffusion.txt", status="old")
64 read(my_unit, settings)
65 close(my_unit)
66 else
67 print *, "No input file reaction_diffusion.txt found; default settings"
68 end if
69
70 call af_add_cc_variable(tree, "u", ix=i_u, n_copies=2)
71 call af_add_cc_variable(tree, "v", ix=i_v, n_copies=2)
72 call af_add_cc_variable(tree, "phi1", ix=i_phi1)
73 call af_add_cc_variable(tree, "phi2", ix=i_phi2)
74 call af_add_cc_variable(tree, "rhs1", ix=i_rhs1)
75 call af_add_cc_variable(tree, "rhs2", ix=i_rhs2)
76 call af_add_cc_variable(tree, "tmp", ix=i_tmp)
77
78 call af_set_cc_methods(tree, i_u, af_bc_neumann_zero)
79 call af_set_cc_methods(tree, i_v, af_bc_neumann_zero)
80
81
82 call af_init(tree, &
83 box_size, &
84 [dtimes(domain_len)], &
85 [dtimes(box_size)], &
86 periodic=periodic)
87
88 call af_print_info(tree)
89
90
91 do
92 call af_loop_box(tree, set_initial_condition)
93 call af_gc_tree(tree, [i_u, i_v])
94 call af_adjust_refinement(tree, &
95 refine_routine, &
96 refine_info)
97 if (refine_info%n_add == 0) exit
98 end do
99
100
101 mg1%i_phi = i_phi1
102 mg1%i_rhs = i_rhs1
103 mg1%i_tmp = i_tmp
104 mg1%sides_bc => af_bc_neumann_zero
105
106 mg2%i_phi = i_phi2
107 mg2%i_rhs = i_rhs2
108 mg2%i_tmp = i_tmp
109 mg2%sides_bc => af_bc_neumann_zero
110
111 output_cnt = 0
112 time = 0
113 it = 0
114 time = 0
115
116 select case (time_integrator)
117 case ("imex")
118 continue
119 case default
120 if (dt > af_min_dr(tree)**2/(2*ndim*max(d1, d2))) then
121 error stop "dt is too large for explicit integration"
122 end if
123 end select
124
125
126 mg1%helmholtz_lambda = 1/(0.5_dp * dt * d1)
127 mg2%helmholtz_lambda = 1/(0.5_dp * dt * d2)
128
129 call mg_init(tree, mg1)
130 call mg_init(tree, mg2)
131
132
133 do
134 it = it + 1
135
136 if (output_cnt * dt_output <= time) then
137 output_cnt = output_cnt + 1
138 write(fname, "(A,I0)") "output/reaction_diffusion_" &
139 // dimname // "_", output_cnt
140 call af_write_silo(tree, trim(fname), output_cnt, time)
141 end if
142
143 select case (time_integrator)
144 case ("imex")
145 call rd_imex(tree)
146 case ("forward_euler")
147 call af_loop_box(tree, forward_euler)
148 case ("midpoint_method")
149 call af_loop_box(tree, midpoint_method_step1)
150 call af_gc_tree(tree, [i_u+1, i_v+1])
151 call af_loop_box(tree, midpoint_method_step2)
152 end select
153
154 call af_gc_tree(tree, [i_u, i_v])
155
156 if (time > end_time) exit
157 time = time + dt
158 end do
159
160contains
161
162 subroutine forward_euler(box)
163 type(box_t), intent(inout) :: box
164 call step_f(box, dt, [i_u, i_v], [i_u, i_v], [i_u, i_v])
165 end subroutine forward_euler
166
167 subroutine midpoint_method_step1(box)
168 type(box_t), intent(inout) :: box
169 call step_f(box, 0.5_dp * dt, [i_u, i_v], [i_u, i_v], [i_u+1, i_v+1])
170 end subroutine midpoint_method_step1
171
172 subroutine midpoint_method_step2(box)
173 type(box_t), intent(inout) :: box
174 call step_f(box, dt, [i_u+1, i_v+1], [i_u, i_v], [i_u, i_v])
175 end subroutine midpoint_method_step2
176
177
178
179
180 subroutine rd_imex(tree)
181 type(af_t), intent(inout) :: tree
182 integer :: n
183 real(dp) :: max_res, max_rhs
184 real(dp), parameter :: max_rel_residual = 1e-6_dp
185
186 call af_loop_box(tree, rd_imex_step1)
187
188 do n = 1, 10
189 call af_tree_maxabs_cc(tree, mg1%i_rhs, max_rhs)
190 call mg_fas_fmg(tree, mg1, .true., .true.)
191
192 call af_tree_maxabs_cc(tree, mg1%i_tmp, max_res)
193 if (max_res/max_rhs < max_rel_residual) exit
194 end do
195
196 do n = 1, 10
197 call af_tree_maxabs_cc(tree, mg2%i_rhs, max_rhs)
198 call mg_fas_fmg(tree, mg2, .true., .true.)
199
200 call af_tree_maxabs_cc(tree, mg2%i_tmp, max_res)
201 if (max_res/max_rhs < max_rel_residual) exit
202 end do
203
204 call af_loop_box(tree, rd_imex_step2)
205 end subroutine rd_imex
206
207 subroutine rd_imex_step1(box)
208 type(box_t), intent(inout) :: box
209 integer :: nc
210
211 nc = box%n_cell
212
213
214 call step_f0(box, dt, [i_u, i_v], [i_u, i_v], [mg1%i_rhs, mg2%i_rhs])
215 call step_f1(box, 0.5_dp * dt, [i_u, i_v], [mg1%i_rhs, mg2%i_rhs], &
216 [mg1%i_rhs, mg2%i_rhs])
217
218 box%cc(dtimes(1:nc), mg1%i_rhs) = -box%cc(dtimes(1:nc), mg1%i_rhs) * &
219 mg1%helmholtz_lambda
220 box%cc(dtimes(1:nc), mg2%i_rhs) = -box%cc(dtimes(1:nc), mg2%i_rhs) * &
221 mg2%helmholtz_lambda
222
223 end subroutine rd_imex_step1
224
225 subroutine rd_imex_step2(box)
226 type(box_t), intent(inout) :: box
227
228 call step_f(box, 0.5_dp * dt, [i_u, i_v], [i_u, i_v], [i_u, i_v])
229 call step_f(box, 0.5_dp * dt, [mg1%i_phi, mg2%i_phi], &
230 [i_u, i_v], [i_u, i_v])
231 end subroutine rd_imex_step2
232
233
234 subroutine step_f0(box, dt, i_deriv, i_prev, i_out)
235 type(box_t), intent(inout) :: box
236 real(dp), intent(in) :: dt
237 integer, intent(in) :: i_deriv(2)
238 integer, intent(in) :: i_prev(2)
239 integer, intent(in) :: i_out(2)
240 real(dp) :: tmp(DTIMES(box%n_cell), 2)
241 integer :: nc
242
243 nc = box%n_cell
244
245
246 select case (equation_type)
247 case ("gs")
248 tmp(dtimes(:), 1) = -box%cc(dtimes(1:nc), i_deriv(1)) * &
249 box%cc(dtimes(1:nc), i_deriv(2))**2 + &
250 gs_f * (1 - box%cc(dtimes(1:nc), i_deriv(1)))
251 tmp(dtimes(:), 2) = box%cc(dtimes(1:nc), i_deriv(1)) * &
252 box%cc(dtimes(1:nc), i_deriv(2))**2 - &
253 (gs_f + gs_k) * box%cc(dtimes(1:nc), i_deriv(2))
254 case ("schnakenberg")
255 tmp(dtimes(:), 1) = kappa * (alpha - box%cc(dtimes(1:nc), i_deriv(1)) + &
256 box%cc(dtimes(1:nc), i_deriv(1))**2 * box%cc(dtimes(1:nc), i_deriv(2)))
257 tmp(dtimes(:), 2) = kappa * (beta - &
258 box%cc(dtimes(1:nc), i_deriv(1))**2 * box%cc(dtimes(1:nc), i_deriv(2)))
259 case default
260 error stop "Invalid equation type"
261 end select
262
263 box%cc(dtimes(1:nc), i_out) = &
264 box%cc(dtimes(1:nc), i_prev) + dt * tmp
265 end subroutine step_f0
266
267
268 subroutine step_f1(box, dt, i_deriv, i_prev, i_out)
269 type(box_t), intent(inout) :: box
270 real(dp), intent(in) :: dt
271 integer, intent(in) :: i_deriv(2)
272 integer, intent(in) :: i_prev(2)
273 integer, intent(in) :: i_out(2)
274 real(dp) :: tmp(DTIMES(box%n_cell), 2)
275 integer :: nc
276
277 nc = box%n_cell
278 call laplacian(box, i_deriv(1), tmp(dtimes(:), 1))
279 call laplacian(box, i_deriv(2), tmp(dtimes(:), 2))
280 tmp(dtimes(:), 1) = tmp(dtimes(:), 1) * d1
281 tmp(dtimes(:), 2) = tmp(dtimes(:), 2) * d2
282
283 box%cc(dtimes(1:nc), i_out) = &
284 box%cc(dtimes(1:nc), i_prev) + dt * tmp
285 end subroutine step_f1
286
287 subroutine step_f(box, dt, i_deriv, i_prev, i_out)
288 type(box_t), intent(inout) :: box
289 real(dp), intent(in) :: dt
290 integer, intent(in) :: i_deriv(2)
291 integer, intent(in) :: i_prev(2)
292 integer, intent(in) :: i_out(2)
293 integer :: nc
294 real(dp) :: tmp(DTIMES(box%n_cell), 2)
295
296 nc = box%n_cell
297
298 call laplacian(box, i_deriv(1), tmp(dtimes(:), 1))
299 call laplacian(box, i_deriv(2), tmp(dtimes(:), 2))
300 tmp(dtimes(:), 1) = tmp(dtimes(:), 1) * d1
301 tmp(dtimes(:), 2) = tmp(dtimes(:), 2) * d2
302
303
304 select case (equation_type)
305 case ("gs")
306 tmp(dtimes(:), 1) = tmp(dtimes(:), 1) - &
307 box%cc(dtimes(1:nc), i_deriv(1)) * &
308 box%cc(dtimes(1:nc), i_deriv(2))**2 + &
309 gs_f * (1 - box%cc(dtimes(1:nc), i_deriv(1)))
310 tmp(dtimes(:), 2) = tmp(dtimes(:), 2) + &
311 box%cc(dtimes(1:nc), i_deriv(1)) * &
312 box%cc(dtimes(1:nc), i_deriv(2))**2 - &
313 (gs_f + gs_k) * box%cc(dtimes(1:nc), i_deriv(2))
314 case ("schnakenberg")
315 tmp(dtimes(:), 1) = tmp(dtimes(:), 1) + &
316 kappa * (alpha - box%cc(dtimes(1:nc), i_deriv(1)) + &
317 box%cc(dtimes(1:nc), i_deriv(1))**2 * box%cc(dtimes(1:nc), i_deriv(2)))
318 tmp(dtimes(:), 2) = tmp(dtimes(:), 2) + &
319 kappa * (beta - &
320 box%cc(dtimes(1:nc), i_deriv(1))**2 * box%cc(dtimes(1:nc), i_deriv(2)))
321 case default
322 error stop "Invalid equation type"
323 end select
324
325 box%cc(dtimes(1:nc), i_out) = box%cc(dtimes(1:nc), i_prev) + dt * tmp
326 end subroutine step_f
327
328
329 subroutine refine_routine(box, cell_flags)
330 type(box_t), intent(in) :: box
331 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
332
333 if (maxval(box%dr) > domain_len/nx_min) then
334 cell_flags = af_do_ref
335 else
336 cell_flags = af_keep_ref
337 end if
338 end subroutine refine_routine
339
340
341 subroutine set_initial_condition(box)
342 type(box_t), intent(inout) :: box
343 integer :: IJK, nc
344 real(dp) :: rr(NDIM), r0(NDIM), u, v
345
346 nc = box%n_cell
347
348 select case (equation_type)
349 case ("gs")
350 do kji_do(0,nc+1)
351 rr = af_r_cc(box, [ijk])
352
353 if (all(abs(rr - 0.5_dp * domain_len) < 10.0_dp/256)) then
354 call random_number(u)
355 call random_number(v)
356 box%cc(ijk, i_u) = 0.5_dp + u_rng_ampl * (u - 0.5_dp)
357 box%cc(ijk, i_v) = 0.25_dp + v_rng_ampl * (v - 0.5_dp)
358 else
359 box%cc(ijk, i_u) = 1.0_dp
360 box%cc(ijk, i_v) = 0.0_dp
361 end if
362 end do; close_do
363 case ("schnakenberg")
364 r0(1) = 1.0_dp/3
365 r0(2:) = 0.5_dp
366
367 do kji_do(0,nc+1)
368 rr = af_r_cc(box, [ijk])
369 call random_number(u)
370 call random_number(v)
371 box%cc(ijk, i_u) = alpha + beta + 1e-3_dp * &
372 exp(-100_dp * sum((rr - r0)**2)) + u_rng_ampl * (u - 0.5_dp)
373 box%cc(ijk, i_v) = beta / (alpha + beta)**2 + v_rng_ampl * (v - 0.5_dp)
374 end do; close_do
375 case default
376 error stop "Unknown equation type"
377 end select
378 end subroutine set_initial_condition
379
380
381 subroutine laplacian(box, i_in, lpl)
382 type(box_t), intent(inout) :: box
383 integer, intent(in) :: i_in
384 real(dp), intent(out) :: lpl(DTIMES(1:box%n_cell))
385 integer :: IJK, nc
386 real(dp) :: idr2(NDIM)
387
388 nc = box%n_cell
389 idr2 = 1 / box%dr**2
390
391 associate(cc => box%cc, n => i_in)
392 do kji_do(1, nc)
393#if NDIM == 1
394 lpl(i) = idr2(1) * (cc(i-1, n) + cc(i+1, n) - 2 * cc(i, n))
395#elif NDIM == 2
396 lpl(i, j) = &
397 idr2(1) * (cc(i-1, j, n) + cc(i+1, j, n) - 2 * cc(i, j, n)) + &
398 idr2(2) * (cc(i, j-1, n) + cc(i, j+1, n) - 2 * cc(i, j, n))
399#elif NDIM == 3
400 lpl(i, j, k) = &
401 idr2(1) * (cc(i-1, j, k, n) + cc(i+1, j, k, n) &
402 - 2 * cc(i, j, k, n)) &
403 + idr2(2) * (cc(i, j-1, k, n) + cc(i, j+1, k, n) &
404 - 2 * cc(i, j, k, n)) &
405 + idr2(3) * (cc(i, j, k-1, n) + cc(i, j, k+1, n) &
406 - 2 * cc(i, j, k, n))
407#endif
408 end do; close_do
409 end associate
410 end subroutine laplacian
411
412end program reaction_diffusion
Module which contains all Afivo modules, so that a user does not have to include them separately.