Commit 8162438a authored by Jannis Teunissen's avatar Jannis Teunissen

Add flag for periodic boundary conditions

parent b2aeb828
......@@ -218,6 +218,9 @@ module m_streamer
! The length of the (square) domain
real(dp), protected :: ST_domain_len = 16e-3_dp
! Whether the domain is periodic in the x (and y)
logical, protected :: ST_periodic = .false.
! Pressure of the gas in bar
real(dp), protected :: ST_gas_pressure = 1.0_dp
......@@ -327,6 +330,8 @@ contains
"The number of grid cells per coordinate in a box")
call CFG_add_get(cfg, "domain_len", ST_domain_len, &
"The length of the domain (m)")
call CFG_add_get(cfg, "periodic", ST_periodic, &
"Whether the domain is periodic")
call CFG_add_get(cfg, "gas_pressure", ST_gas_pressure, &
"The gas pressure (bar), used for photoionization")
call CFG_add_get(cfg, "gas_frac_O2", ST_gas_frac_O2, &
......
......@@ -278,6 +278,7 @@ contains
real(dp) :: dr
integer :: id
integer :: ix_list($D, 1) ! Spatial indices of initial boxes
integer :: nb_list(2*$D, 1) ! Index of neighbors
integer :: n_boxes_init = 1000
dr = ST_domain_len / ST_box_size
......@@ -296,8 +297,20 @@ contains
id = 1 ! One box ...
ix_list(:, id) = 1 ! With index 1,1 ...
! Create the base mesh
call a$D_set_base(tree, 1, ix_list)
if (ST_periodic) then
if (ST_cylindrical) error stop "Cannot have periodic cylindrical domain"
nb_list(:, 1) = -1
nb_list(a$D_neighb_lowx, 1) = 1
nb_list(a$D_neighb_highx, 1) = 1
#if $D == 3
nb_list(a$D_neighb_lowy, 1) = 1
nb_list(a$D_neighb_highy, 1) = 1
#endif
call a$D_set_base(tree, 1, ix_list, nb_list)
else
! Create the base mesh
call a$D_set_base(tree, 1, ix_list)
end if
end subroutine init_tree
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment