A simplified model for ionization waves and/or streamers in 2D.
1
2
3
4program simple_streamer
5
15
16 implicit none
17
18 integer :: i, n
19 character(len=100) :: fname
20 type(af_t) :: tree
21 type(mg_t) :: mg
22 type(ref_info_t) :: refine_info
23
24
25 integer :: i_elec
26 integer :: i_pion
27 integer :: i_elec_old
28 integer :: i_pion_old
29 integer :: i_phi
30 integer :: i_fld
31 integer :: i_rhs
32
33
34 integer :: f_elec
35 integer :: f_fld
36
37
38 real(dp), parameter :: end_time = 3e-9_dp
39 real(dp), parameter :: dt_output = 20e-11_dp
40 real(dp), parameter :: dt_max = 20e-11_dp
41 integer, parameter :: ref_per_steps = 2
42 integer, parameter :: box_size = 8
43
44
45 real(dp), parameter :: applied_field = -0.8e7_dp
46 real(dp), parameter :: mobility = 0.03_dp
47 real(dp), parameter :: diffusion_c = 0.2_dp
48
49
50 real(dp), parameter :: domain_length = 2e-3_dp
51 real(dp), parameter :: refine_max_dx = 1e-3_dp
52 real(dp), parameter :: refine_min_dx = 1e-9_dp
53
54
55 real(dp), parameter :: init_density = 1e15_dp
56 real(dp), parameter :: init_y_min = 0.8_dp * domain_length
57 real(dp), parameter :: init_y_max = 0.9_dp * domain_length
58
59
60 real(dp) :: dt
61 real(dp) :: time
62 integer :: output_count
63
64
65 real(dp) :: sum_elec, sum_pion
66
67 call af_add_cc_variable(tree, "elec", ix=i_elec, n_copies=2)
68 i_elec_old = af_find_cc_variable(tree, "elec_2")
69 call af_add_cc_variable(tree, "pion", ix=i_pion, n_copies=2)
70 i_pion_old = af_find_cc_variable(tree, "pion_2")
71 call af_add_cc_variable(tree, "phi", ix=i_phi)
72 call af_add_cc_variable(tree, "fld", ix=i_fld)
73 call af_add_cc_variable(tree, "rhs", ix=i_rhs)
74
75 call af_add_fc_variable(tree, "f_elec", ix=f_elec)
76 call af_add_fc_variable(tree, "f_fld", ix=f_fld)
77
78
79 call af_init(tree, &
80 box_size, &
81 [domain_length, domain_length], &
82 [box_size, box_size], &
83 periodic=[.true., .false.])
84
85
86
87 mg%i_phi = i_phi
88 mg%i_tmp = i_fld
89 mg%i_rhs = i_rhs
90
91
92 mg%sides_bc => sides_bc_pot
93
94
95 call mg_init(tree, mg)
96
97 call af_set_cc_methods(tree, i_elec, af_bc_dirichlet_zero, &
98 prolong=af_prolong_limit)
99 call af_set_cc_methods(tree, i_pion, af_bc_dirichlet_zero, &
100 prolong=af_prolong_limit)
101 call af_set_cc_methods(tree, i_fld, af_bc_neumann_zero)
102
103 output_count = 0
104 time = 0
105
106
107 do
108
109 call af_loop_box(tree, set_initial_condition)
110
111
112
113
114 call compute_fld(tree, .false.)
115
116
117
118
119
120
121 call af_adjust_refinement(tree, &
122 refinement_routine, &
123 refine_info)
124
125
126 if (refine_info%n_add == 0 .and. refine_info%n_rm == 0) exit
127 end do
128
129 call af_print_info(tree)
130
131 do
132
133 call af_reduction(tree, &
134 max_dt, &
135 get_min, &
136 dt_max, &
137 dt)
138
139 if (dt < 1e-14) then
140 print *, "dt getting too small, instability?", dt
141 time = end_time + 1.0_dp
142 end if
143
144
145 if (output_count * dt_output <= time) then
146 output_count = output_count + 1
147 write(fname, "(A,I6.6)") "output/simple_streamer_", output_count
148
149
150
151 call af_write_silo(tree, fname, output_count, time)
152
153 call af_tree_sum_cc(tree, i_elec, sum_elec)
154 call af_tree_sum_cc(tree, i_pion, sum_pion)
155 print *, "sum(pion-elec)/sum(pion)", (sum_pion - sum_elec)/sum_pion
156 end if
157
158 if (time > end_time) exit
159
160
161 do n = 1, ref_per_steps
162 time = time + dt
163
164
165 call af_tree_copy_cc(tree, i_elec, i_elec_old)
166 call af_tree_copy_cc(tree, i_pion, i_pion_old)
167
168
169 do i = 1, 2
170
171 call af_loop_tree(tree, fluxes_koren, leaves_only=.true.)
172 call af_consistent_fluxes(tree, [f_elec])
173
174
175 call af_loop_box_arg(tree, update_solution, [dt], &
176 leaves_only=.true.)
177
178
179 if (i == 1) call compute_fld(tree, .true.)
180 end do
181
182
183 call af_loop_box(tree, average_dens)
184
185
186 call compute_fld(tree, .true.)
187 end do
188
189
190 call af_gc_tree(tree, [i_elec, i_pion])
191
192
193
194
195
196
197
198 call af_adjust_refinement(tree, &
199 refinement_routine, &
200 refine_info, &
201 4)
202
203 if (refine_info%n_add > 0 .or. refine_info%n_rm > 0) then
204
205 call compute_fld(tree, .true.)
206 end if
207 end do
208
209contains
210
211
212 subroutine refinement_routine(box, cell_flags)
213 type(box_t), intent(in) :: box
214 integer, intent(out) :: cell_flags(box%n_cell, box%n_cell)
215 integer :: i, j, nc
216 real(dp) :: dx, dens, fld, adx, xy(2)
217
218 nc = box%n_cell
219 dx = maxval(box%dr)
220
221 do j = 1, nc
222 do i = 1, nc
223 xy = af_r_cc(box, [i,j])
224 dens = box%cc(i, j, i_elec)
225 fld = box%cc(i, j, i_fld)
226 adx = get_alpha(fld) * dx
227
228 if (dens > 1e0_dp .and. adx > 0.8_dp) then
229 cell_flags(i, j) = af_do_ref
230 else if (dx < 1.25e-5_dp .and. adx < 0.1_dp) then
231 cell_flags(i, j) = af_rm_ref
232 else
233 cell_flags(i, j) = af_keep_ref
234 end if
235 end do
236 end do
237
238
239 if (dx > refine_max_dx) then
240 cell_flags = af_do_ref
241 else if (dx < 2 * refine_min_dx) then
242 where(cell_flags == af_do_ref) cell_flags = af_keep_ref
243 else if (dx > 0.5_dp * refine_max_dx) then
244 where(cell_flags == af_rm_ref) cell_flags = af_keep_ref
245 end if
246
247 end subroutine refinement_routine
248
249
250 subroutine set_initial_condition(box)
251 type(box_t), intent(inout) :: box
252 integer :: i, j, nc
253 real(dp) :: xy(2), normal_rands(2), vol
254
255 nc = box%n_cell
256
257 do j = 0, nc+1
258 do i = 0, nc+1
259 xy = af_r_cc(box, [i,j])
260 vol = (box%dr(1) * box%dr(2))**1.5_dp
261
262 if (xy(2) > init_y_min .and. xy(2) < init_y_max) then
263
264 normal_rands = two_normals(vol * init_density, &
265 sqrt(vol * init_density))
266
267 box%cc(i, j, i_elec) = abs(normal_rands(1)) / vol
268 else
269 box%cc(i, j, i_elec) = 0
270 end if
271 end do
272 end do
273
274 box%cc(:, :, i_pion) = box%cc(:, :, i_elec)
275 box%cc(:, :, i_phi) = 0
276
277 end subroutine set_initial_condition
278
279
280
281 function two_normals(mean, sigma) result(rands)
282 real(dp), intent(in) :: mean, sigma
283 real(dp) :: rands(2), sum_sq
284
285 do
286 call random_number(rands)
287 rands = 2 * rands - 1
288 sum_sq = sum(rands**2)
289 if (sum_sq < 1.0_dp .and. sum_sq > 0.0_dp) exit
290 end do
291 rands = rands * sqrt(-2 * log(sum_sq) / sum_sq)
292 rands = mean + rands * sigma
293 end function two_normals
294
295
296 real(dp) function get_min(a, b)
297 real(dp), intent(in) :: a, b
298 get_min = min(a, b)
299 end function get_min
300
301
302 real(dp) function max_dt(box)
303 type(box_t), intent(in) :: box
304 real(dp), parameter :: UC_eps0 = 8.8541878176d-12
305 real(dp), parameter :: UC_elem_charge = 1.6022d-19
306 integer :: i, j, nc
307 real(dp) :: fld(2)
308 real(dp) :: dt_cfl, dt_dif, dt_drt
309
310 nc = box%n_cell
311 dt_cfl = dt_max
312
313 do j = 1, nc
314 do i = 1, nc
315 fld(1) = 0.5_dp * (box%fc(i, j, 1, f_fld) + &
316 box%fc(i+1, j, 1, f_fld))
317 fld(2) = 0.5_dp * (box%fc(i, j, 2, f_fld) + &
318 box%fc(i, j+1, 2, f_fld))
319
320
321 dt_cfl = min(dt_cfl, 0.5_dp / sum(abs(fld * mobility) / box%dr))
322 end do
323 end do
324
325
326 dt_drt = uc_eps0 / (uc_elem_charge * mobility * &
327 maxval(box%cc(1:nc, 1:nc, i_elec)) + epsilon(1.0_dp))
328
329
330 dt_dif = 0.25_dp * minval(box%dr)**2 / max(diffusion_c, epsilon(1.0_dp))
331
332 max_dt = min(dt_cfl, dt_drt, dt_dif)
333 end function max_dt
334
335
336
337
338
339
340
341 elemental function get_alpha(fld) result(alpha)
342 real(dp), intent(in) :: fld
343 real(dp) :: alpha
344 real(dp), parameter :: c0 = 1.04e1_dp
345 real(dp), parameter :: c1 = 0.601_dp
346 real(dp), parameter :: c2 = 1.86e7_dp
347
348 alpha = exp(c0) * (abs(fld) * 1e-5_dp)**c1 * exp(-c2 / abs(fld))
349 end function get_alpha
350
351
352
353 subroutine compute_fld(tree, have_guess)
354 type(af_t), intent(inout) :: tree
355 logical, intent(in) :: have_guess
356 real(dp), parameter :: fac = 1.6021766208e-19_dp / 8.8541878176e-12_dp
357 integer :: lvl, i, id, nc
358
359 nc = tree%n_cell
360
361
362
363 do lvl = 1, tree%highest_lvl
364
365 do i = 1, size(tree%lvls(lvl)%leaves)
366 id = tree%lvls(lvl)%leaves(i)
367 tree%boxes(id)%cc(:, :, i_rhs) = fac * (&
368 tree%boxes(id)%cc(:, :, i_elec) - &
369 tree%boxes(id)%cc(:, :, i_pion))
370 end do
371
372 end do
373
374
375
376 call mg_fas_fmg(tree, mg, .false., have_guess)
377
378
379 call af_loop_box(tree, fld_from_pot)
380
381
382 call af_gc_tree(tree, [i_fld])
383 end subroutine compute_fld
384
385
386 subroutine fld_from_pot(box)
387 type(box_t), intent(inout) :: box
388 integer :: nc
389 real(dp) :: inv_dr(2)
390
391 nc = box%n_cell
392 inv_dr = 1 / box%dr
393
394 box%fc(1:nc+1, 1:nc, 1, f_fld) = inv_dr(1) * &
395 (box%cc(0:nc, 1:nc, i_phi) - box%cc(1:nc+1, 1:nc, i_phi))
396 box%fc(1:nc, 1:nc+1, 2, f_fld) = inv_dr(2) * &
397 (box%cc(1:nc, 0:nc, i_phi) - box%cc(1:nc, 1:nc+1, i_phi))
398
399 box%cc(1:nc, 1:nc, i_fld) = sqrt(&
400 0.25_dp * (box%fc(1:nc, 1:nc, 1, f_fld) + &
401 box%fc(2:nc+1, 1:nc, 1, f_fld))**2 + &
402 0.25_dp * (box%fc(1:nc, 1:nc, 2, f_fld) + &
403 box%fc(1:nc, 2:nc+1, 2, f_fld))**2)
404 end subroutine fld_from_pot
405
406
407 subroutine fluxes_koren(tree, id)
409 type(af_t), intent(inout) :: tree
410 integer, intent(in) :: id
411 real(dp) :: inv_dr(2)
412 real(dp) :: cc(-1:tree%n_cell+2, -1:tree%n_cell+2, 1)
413 real(dp), allocatable :: v(:, :, :), dc(:, :, :)
414 integer :: nc
415
416 nc = tree%n_cell
417 inv_dr = 1/tree%boxes(id)%dr
418
419 allocate(v(1:nc+1, 1:nc+1, 2))
420 allocate(dc(1:nc+1, 1:nc+1, 2))
421
422 call af_gc2_box(tree, id, [i_elec], cc)
423
424 v = -mobility * tree%boxes(id)%fc(:, :, :, f_fld)
425 dc = diffusion_c
426
427 call flux_koren_2d(cc(:, :, 1), v, nc, 2)
428 call flux_diff_2d(cc(:, :, 1), dc, inv_dr, nc, 2)
429
430 tree%boxes(id)%fc(:, :, :, f_elec) = v + dc
431 end subroutine fluxes_koren
432
433
434 subroutine average_dens(box)
435 type(box_t), intent(inout) :: box
436 box%cc(:, :, i_elec) = 0.5_dp * (box%cc(:, :, i_elec) + box%cc(:, :, i_elec_old))
437 box%cc(:, :, i_pion) = 0.5_dp * (box%cc(:, :, i_pion) + box%cc(:, :, i_pion_old))
438 end subroutine average_dens
439
440
441 subroutine update_solution(box, dt)
442 type(box_t), intent(inout) :: box
443 real(dp), intent(in) :: dt(:)
444 real(dp) :: inv_dr(2), src, sflux, fld
445 real(dp) :: alpha
446 integer :: i, j, nc
447
448 nc = box%n_cell
449 inv_dr = 1/box%dr
450
451 do j = 1, nc
452 do i = 1, nc
453 fld = box%cc(i,j, i_fld)
454 alpha = get_alpha(fld)
455 src = box%cc(i, j, i_elec) * mobility * abs(fld) * alpha
456
457 sflux = inv_dr(1) * (box%fc(i, j, 1, f_elec) - &
458 box%fc(i+1, j, 1, f_elec)) + &
459 inv_dr(2) * (box%fc(i, j, 2, f_elec) - &
460 box%fc(i, j+1, 2, f_elec))
461
462 box%cc(i, j, i_elec) = box%cc(i, j, i_elec) + (src + sflux) * dt(1)
463 box%cc(i, j, i_pion) = box%cc(i, j, i_pion) + src * dt(1)
464 end do
465 end do
466
467 end subroutine update_solution
468
469
470 subroutine sides_bc_pot(box, nb, iv, coords, bc_val, bc_type)
471 type(box_t), intent(in) :: box
472 integer, intent(in) :: nb
473 integer, intent(in) :: iv
474 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
475 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
476 integer, intent(out) :: bc_type
477
478 select case (nb)
479 case (af_neighb_lowy)
480 bc_type = af_bc_neumann
481 bc_val = applied_field
482 case (af_neighb_highy)
483 bc_type = af_bc_dirichlet
484 bc_val = 0.0_dp
485 case default
486 stop "sides_bc_pot: unspecified boundary"
487 end select
488 end subroutine sides_bc_pot
489
490end program simple_streamer
This module contains the core routines of Afivo, namely those that deal with initializing and changin...
Module containing a couple flux schemes for solving hyperbolic problems explicitly,...
This module contains routines related to the filling of ghost cells. Note that corner ghost cells are...
This module contains the geometric multigrid routines that come with Afivo.
This module contains routines for writing output files with Afivo. The Silo format should probably be...
This module contains the routines related to prolongation: going from coarse to fine variables.
This module contains routines for restriction: going from fine to coarse variables.
This module contains the basic types and constants that are used in the NDIM-dimensional version of A...
This module contains all kinds of different 'helper' routines for Afivo. If the number of routines fo...
This module contains wrapper functions to simplify writing Silo files.