Afivo 0.3
All Classes Namespaces Functions Variables Pages
simple_streamer.f90

A simplified model for ionization waves and/or streamers in 2D.

A simplified model for ionization waves and/or streamers in 2D

1!> \example simple_streamer.f90
2!>
3!> A simplified model for ionization waves and/or streamers in 2D
4program simple_streamer
5
6 use m_af_types
7 use m_af_core
9 use m_af_utils
12 use m_af_output
13 use m_write_silo
14 use m_af_prolong
15
16 implicit none
17
18 integer :: i, n
19 character(len=100) :: fname
20 type(af_t) :: tree ! This contains the full grid information
21 type(mg_t) :: mg ! Multigrid option struct
22 type(ref_info_t) :: refine_info
23
24 ! Indices of cell-centered variables
25 integer :: i_elec ! Electron density
26 integer :: i_pion ! Positive ion density
27 integer :: i_elec_old ! For time-stepping scheme
28 integer :: i_pion_old ! For time-stepping scheme
29 integer :: i_phi ! Electrical potential
30 integer :: i_fld ! Electric field norm
31 integer :: i_rhs ! Source term Poisson
32
33 ! Indices of face-centered variables **
34 integer :: f_elec ! Electron flux
35 integer :: f_fld ! Electric field vector
36
37 ! Simulation parameters
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 ! Physical parameters
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 ! Computational domain
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 ! Settings for the initial conditions
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 ! Simulation variables
60 real(dp) :: dt
61 real(dp) :: time
62 integer :: output_count
63
64 ! To test charge conservation
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 ! Initialize the tree (which contains all the mesh information)
79 call af_init(tree, & ! Tree to initialize
80 box_size, & ! A box contains box_size**DIM cells
81 [domain_length, domain_length], &
82 [box_size, box_size], &
83 periodic=[.true., .false.])
84
85
86 ! Set the multigrid options. First define the variables to use
87 mg%i_phi = i_phi
88 mg%i_tmp = i_fld
89 mg%i_rhs = i_rhs
90
91 ! Routines to use for ...
92 mg%sides_bc => sides_bc_pot ! Filling ghost cell on physical boundaries
93
94 ! This routine always needs to be called when using multigrid
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 ! Number of output files written
104 time = 0 ! Simulation time (all times are in s)
105
106 ! Set up the initial conditions
107 do
108 ! For each box in tree, set the initial conditions
109 call af_loop_box(tree, set_initial_condition)
110
111 ! Compute electric field on the tree.
112 ! First perform multigrid to get electric potential,
113 ! then take numerical gradient to geld field.
114 call compute_fld(tree, .false.)
115
116 ! Adjust the refinement of a tree using refine_routine (see below) for grid
117 ! refinement.
118 ! Routine af_adjust_refinement sets the bit af_bit_new_children for each box
119 ! that is refined. On input, the tree should be balanced. On output,
120 ! the tree is still balanced, and its refinement is updated (with at most
121 call af_adjust_refinement(tree, & ! tree
122 refinement_routine, & ! Refinement function
123 refine_info) ! Information about refinement
124
125 ! If no new boxes have been added or removed, exit the loop
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 ! Get a new time step, which is at most dt_max
133 call af_reduction(tree, & ! Tree to do the reduction on
134 max_dt, & ! function
135 get_min, & ! function
136 dt_max, & ! Initial value for the reduction
137 dt) ! Result of the reduction
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 ! Every dt_output, write output
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 ! Write the cell centered data of a tree to a Silo file. Only the
150 ! leaves of the tree are used
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 ! We perform n_steps between mesh-refinements
161 do n = 1, ref_per_steps
162 time = time + dt
163
164 ! Copy previous solution
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 ! Two forward Euler steps over dt
169 do i = 1, 2
170 ! First calculate fluxes
171 call af_loop_tree(tree, fluxes_koren, leaves_only=.true.)
172 call af_consistent_fluxes(tree, [f_elec])
173
174 ! Update the solution
175 call af_loop_box_arg(tree, update_solution, [dt], &
176 leaves_only=.true.)
177
178 ! Compute new field on first iteration
179 if (i == 1) call compute_fld(tree, .true.)
180 end do
181
182 ! Take average of phi_old and phi (explicit trapezoidal rule)
183 call af_loop_box(tree, average_dens)
184
185 ! Compute field with new density
186 call compute_fld(tree, .true.)
187 end do
188
189 ! Fill ghost cells for i_pion and i_elec
190 call af_gc_tree(tree, [i_elec, i_pion])
191
192 ! Adjust the refinement of a tree using refine_routine (see below) for grid
193 ! refinement.
194 ! Routine af_adjust_refinement sets the bit af_bit_new_children for each box
195 ! that is refined. On input, the tree should be balanced. On output,
196 ! the tree is still balanced, and its refinement is updated (with at most
197 ! one level per call).
198 call af_adjust_refinement(tree, & ! tree
199 refinement_routine, & ! Refinement function
200 refine_info, & ! Information about refinement
201 4) ! Buffer width (in cells)
202
203 if (refine_info%n_add > 0 .or. refine_info%n_rm > 0) then
204 ! Compute the field on the new mesh
205 call compute_fld(tree, .true.)
206 end if
207 end do
208
209contains
210
211 !> This routine sets the refinement flag for boxes(id)
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 ! Make sure we don't have or get a too fine or too coarse grid
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 !> This routine sets the initial conditions for each box
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 ! Approximate Poisson distribution with normal distribution
264 normal_rands = two_normals(vol * init_density, &
265 sqrt(vol * init_density))
266 ! Prevent negative numbers
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 ! Inital potential set to zero
276
277 end subroutine set_initial_condition
278
279 !> Return two normal random variates
280 !> http://en.wikipedia.org/wiki/Marsaglia_polar_method
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 !> This function computes the minimum val of a and b
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 !> Get maximum time step based on e.g. CFL criteria
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 ! The 0.5 is here because of the explicit trapezoidal rule
321 dt_cfl = min(dt_cfl, 0.5_dp / sum(abs(fld * mobility) / box%dr))
322 end do
323 end do
324
325 ! Dielectric relaxation time
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 ! Diffusion condition
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 !> This function gets the alpha value
337 !>
338 !> Taken from: Spatially hybrid computations for streamer discharges: II. Fully
339 !> 3D simulations, Chao Li, Ute Ebert, Willem Hundsdorfer, J. Comput. Phys.
340 !> 231, 1020-1050 (2012), doi:10.1016/j.jcp.2011.07.023
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 ! Compute electric field on the tree. First perform multigrid to get electric
352 ! potential, then take numerical gradient to geld field.
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 ! Set the source term (rhs)
362 !$omp parallel private(lvl, i, id)
363 do lvl = 1, tree%highest_lvl
364 !$omp do
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 !$omp end do nowait
372 end do
373 !$omp end parallel
374
375 ! Perform an FMG cycle (logicals: store residual, first call)
376 call mg_fas_fmg(tree, mg, .false., have_guess)
377
378 ! Compute field from potential
379 call af_loop_box(tree, fld_from_pot)
380
381 ! Set the field norm also in ghost cells
382 call af_gc_tree(tree, [i_fld])
383 end subroutine compute_fld
384
385 ! Compute electric field from electrical potential
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 ! Compute the electron fluxes due to drift and diffusion
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 ! Take average of new and old electron/ion density for explicit trapezoidal rule
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 ! Advance solution over dt based on the fluxes / source term, using forward Euler
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 !> This fills ghost cells near physical boundaries for the potential
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...
Definition: m_af_core.f90:3
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...
Definition: m_af_output.f90:3
This module contains the routines related to prolongation: going from coarse to fine variables.
Definition: m_af_prolong.f90:3
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...
Definition: m_af_types.f90:3
This module contains all kinds of different 'helper' routines for Afivo. If the number of routines fo...
Definition: m_af_utils.f90:4
This module contains wrapper functions to simplify writing Silo files.
Definition: m_write_silo.f90:3