m_streamer.f90

Go to the documentation of this file.

!> This module contains several pre-defined variables like:
!! * Indices of cell-centered variables
!! * Names of the cell-centered variables
!! * Indices of face-centered variables
!! * Indices of transport data
module m_streamer
#include "../afivo/src/cpp_macros.h"
  use m_types
  use m_af_all
  use m_random
  use m_lookup_table
  use m_config
 
  implicit none
  private
 
  !> Index of electrical potential
  integer, public, protected :: i_phi          = -1
  !> Index of electron density
  integer, public, protected :: i_electron     = -1
  !> Index of electron density (in species list)
  integer, public, protected :: ix_electron    = -1
  !> Index of electron energy density
  integer, public, protected :: i_electron_energy = -1
  !> Index of first positive ion species
  integer, public, protected :: i_1pos_ion     = -1
  !> Index of first positive ion (in species list)
  integer, public, protected :: ix_1pos_ion    = -1
  !> Index of electric field norm
  integer, public, protected :: i_electric_fld = -1
  !> Index of source term Poisson
  integer, public, protected :: i_rhs          = -1
  !> Index of temporary variable
  integer, public, protected :: i_tmp          = -1
  !> Index can be set to include a dielectric
  integer, public, protected :: i_eps          = -1
  !> Index can be set to include an electrode
  integer, public, protected :: i_lsf          = -1
 
  !> Index of all densities that evolve in time
  integer, public, protected, allocatable :: all_densities(:)
 
  !> Include deposited power density in output
  logical, public, protected :: compute_power_density = .false.
  !> Index of deposited power density
  integer, public, protected :: i_power_density = -1
 
  !> Index of correction factor for source terms
  integer, public, protected :: i_srcfac = -1
 
  !> Index of electron flux
  integer, public, protected :: flux_elec    = -1
  !> Index of electron energy flux
  integer, public, protected :: flux_energy  = -1
  !> Index of electric field vector
  integer, public, protected :: electric_fld = -1
 
  !> Number of flux variables
  integer, public, protected :: flux_num_species = -1
  !> Number of electron flux variables
  integer, public, protected :: flux_num_electron_vars = -1
  !> List of all flux variables (face-centered index)
  integer, public, protected, allocatable :: flux_variables(:)
  !> List of all flux species (cell-centered index)
  integer, public, protected, allocatable :: flux_species(:)
  !> List of the charges of the flux species
  integer, public, protected, allocatable :: flux_species_charge(:)
  !> List of the signs of the charges of the flux species (+- 1)
  integer, public, protected, allocatable :: flux_species_charge_sign(:)
  !> List of positive ion fluxes (useful for secondary emission)
  integer, public, protected, allocatable :: flux_pos_ion(:)
 
  !> Whether cylindrical coordinates are used
  logical, public, protected :: st_cylindrical = .false.
 
  !> Whether a dielectric is used
  logical, public, protected :: st_use_dielectric = .false.
 
  !> Whether to include an electrode
  logical, public, protected :: st_use_electrode = .false.
 
  !> Boundary condition for the plasma species
  procedure(af_subr_bc), public, protected, pointer :: &
       bc_species => null()
 
  !> Multigrid option structure
  type(mg_t), public :: mg
 
  !> Random number generator
  type(rng_t), public :: st_rng
 
  !> Parallel random number generator
  type(prng_t), public :: st_prng
 
  !> Avoid dielectric relaxation time step constraint by limiting flux
  logical, public, protected :: st_drt_limit_flux = .false.
 
  !> Ensure that flux limiting does not lead to fields higher than this
  real(dp), public, protected :: st_drt_max_field = 1.0e100_dp
 
  !> Use source factor to prevent unphysical effects due to diffusion
  integer, public, protected :: st_source_factor
 
  integer, public, parameter :: source_factor_none = 0
  integer, public, parameter :: source_factor_flux = 1
  integer, public, parameter :: source_factor_original_flux = 2
 
  !> Minimum number of electrons per cell to include source terms
  real(dp), public, protected :: st_source_min_electrons_per_cell = -1e100_dp
 
  !> End time of the simulation
  real(dp), public, protected :: st_end_time = 10e-9_dp
 
  !> Whether streamer length is used as a simulation stopping
  logical, public, protected :: st_use_end_streamer_length = .false.
 
  !> Wait n steps before initializing streamer begin position
  integer, public, protected :: st_initial_streamer_pos_steps_wait = 5
 
  !> Streamer length at which the simulation will stop
  real(dp), public, protected :: st_end_streamer_length = 15e-3
 
  !> The size of the boxes that we use to construct our mesh
  integer, public, protected :: st_box_size = 8
 
  !> Size of the coarse grid
  integer, public, protected :: st_coarse_grid_size(ndim) = -1
 
  !> Domain length per dimension
  real(dp), public, protected :: st_domain_len(ndim) = 16e-3_dp
 
  !> Origin of domain
  real(dp), public, protected :: st_domain_origin(ndim) = 0.0_dp
 
  !> Whether the domain is periodic (per dimension)
  logical, public, protected :: st_periodic(ndim) = .false.
 
  !> Whether to limit plasma reactions to a certain region
  logical, public, protected :: st_plasma_region_enabled = .false.
 
  !> Limit plasma reactions to coordinates between rmin and rmax
  real(dp), public, protected :: st_plasma_region_rmin(ndim) = -1e100_dp
 
  !> Limit plasma reactions to coordinates between rmin and rmax
  real(dp), public, protected :: st_plasma_region_rmax(ndim) = 1e100_dp
 
  !> Number of V-cycles to perform per time step
  integer, public, protected :: st_multigrid_num_vcycles = 2
 
  ! Stop multigrid when residual is smaller than this factor times max(|rhs|)
  real(dp), public, protected :: st_multigrid_max_rel_residual = 1e-4_dp
 
  !> Global time
  real(dp), public :: global_time = 0.0_dp
 
  !> Global time step
  real(dp), public :: global_dt = 0.0_dp
 
  !> Current sum of reaction rates per thread
  real(dp), public, allocatable :: st_current_rates(:, :)
 
  !> Global sum of reaction rates
  real(dp), public, allocatable :: st_global_rates(:)
 
  !> Current sum of J.E per thread
  real(dp), public, allocatable :: st_current_jdote(:, :)
 
  !> Per how many iterations the electric current is computed
  integer, public, protected :: current_update_per_steps = 1000*1000
 
  !> Electric current through electrodes due to J.E
  real(dp), public :: st_global_jdote_current
 
  !> Electric current through electrodes due to displacement current
  real(dp), public :: st_global_displ_current
 
  !> Global sum of J.E
  real(dp), public :: st_global_jdote
 
  ! To keep track of the computational cost of different parts
  real(dp), public :: wc_time_flux = 0.0_dp
  real(dp), public :: wc_time_source = 0.0_dp
  real(dp), public :: wc_time_copy_state = 0.0_dp
  real(dp), public :: wc_time_field = 0.0_dp
  real(dp), public :: wc_time_output = 0.0_dp
  real(dp), public :: wc_time_refine = 0.0_dp
  real(dp), public :: wc_time_photoi = 0.0_dp
 
  !> Method used to prolong (interpolate) densities
  procedure(af_subr_prolong), pointer, public, protected :: &
       st_prolongation_method => null()
 
  public :: st_initialize
 
contains
 
  !> Create the configuration file with default values
  subroutine st_initialize(tree, cfg, ndim)
    use iso_fortran_env, only: int64
    use m_config
    use omp_lib
    use m_chemistry
    use m_units_constants
    use m_gas
    use m_transport_data
    use m_dt
    use m_model
    type(af_t), intent(inout)  :: tree
    type(cfg_t), intent(inout) :: cfg  !< The configuration for the simulation
    integer, intent(in)        :: ndim !< Number of dimensions
    integer                    :: n, k, n_threads, ix_chemistry
    character(len=name_len)    :: prolong_method, bc_method
    character(len=name_len)    :: source_factor = "none"
    character(len=string_len)  :: tmp_str
    integer                    :: rng_int4_seed(4) = &
         [8123, 91234, 12399, 293434]
    integer(int64)             :: rng_int8_seed(2)
    real(dp)                   :: tmp
    integer                    :: flux_ix
    logical                    :: write_source_factor = .false.
 
    ! Set index of electrons
    i_electron = af_find_cc_variable(tree, "e")
    ix_electron = species_index("e")
 
    ! Set index of first positive ion species
    do n = n_gas_species+1, n_species
       if (species_charge(n) == 1) then
          i_1pos_ion = species_itree(n)
          ix_1pos_ion = n
          exit
       end if
    end do
 
    if (i_1pos_ion == -1) error stop "No positive ion species (1+) found"
 
    ! Set flux species
    call af_add_fc_variable(tree, "flux_elec", ix=flux_elec, &
         write_binary=.false.)
    call af_add_fc_variable(tree, "field", ix=electric_fld)
 
    all_densities = species_itree(n_gas_species+1:n_species)
 
    if (model_has_energy_equation) then
       i_electron_energy = af_find_cc_variable(tree, "e_energy")
       call af_add_fc_variable(tree, "flux_energy", ix=flux_energy, &
            write_binary=.false.)
       flux_num_electron_vars = 2
    else
       flux_num_electron_vars = 1
    end if
 
    flux_num_species = flux_num_electron_vars+transport_data_ions%n_mobile_ions
    allocate(flux_species(flux_num_species))
    allocate(flux_species_charge(flux_num_species))
    allocate(flux_species_charge_sign(flux_num_species))
    allocate(flux_variables(flux_num_species))
 
    flux_species(1)             = i_electron
    flux_species_charge(1)      = -1
    flux_species_charge_sign(1) = -1
    flux_variables(1)           = flux_elec
 
    if (model_has_energy_equation) then
       flux_species(2) = i_electron_energy
       flux_species_charge(2) = 0
       flux_species_charge_sign(2) = -1 ! Used to determine upwind direction
       flux_variables(2) = flux_energy
    end if
 
    do n = 1, transport_data_ions%n_mobile_ions
       flux_ix = flux_num_electron_vars + n
       flux_species(flux_ix) = af_find_cc_variable(tree, &
            trim(transport_data_ions%names(n)))
 
       ! Get index in chemistry list and determine charge
       ix_chemistry = species_index(trim(transport_data_ions%names(n)))
       flux_species_charge(flux_ix) = species_charge(ix_chemistry)
       flux_species_charge_sign(flux_ix) = sign(1, species_charge(ix_chemistry))
 
       call af_add_fc_variable(tree, trim(transport_data_ions%names(n)), &
            ix=flux_variables(flux_ix), write_binary=.false.)
    end do
 
    ! Create a list of positive ion fluxes for secondary emission
    n = count(flux_species_charge > 0)
    allocate(flux_pos_ion(n))
 
    k = 0
    do n = 1, size(flux_species_charge)
       if (flux_species_charge(n) > 0) then
          k = k + 1
          flux_pos_ion(k) = flux_variables(n)
       end if
    end do
 
    ! Add one copy so that the old value can be restored
    call af_add_cc_variable(tree, "phi", ix=i_phi, n_copies=2)
    call af_add_cc_variable(tree, "electric_fld", ix=i_electric_fld)
    call af_add_cc_variable(tree, "rhs", ix=i_rhs)
    call af_add_cc_variable(tree, "tmp", write_out=.false., &
         write_binary=.false., ix=i_tmp)
 
    call cfg_add_get(cfg, "cylindrical", st_cylindrical, &
         "Whether cylindrical coordinates are used (only in 2D)")
 
    call cfg_add_get(cfg, "use_dielectric", st_use_dielectric, &
         "Whether a dielectric is used (experimental)")
    if (st_use_dielectric) then
       call af_add_cc_variable(tree, "eps", ix=i_eps)
       call af_set_cc_methods(tree, i_eps, af_bc_neumann_zero, &
            af_gc_prolong_copy, af_prolong_zeroth)
    end if
 
    call cfg_add_get(cfg, "use_electrode", st_use_electrode, &
         "Whether to include an electrode")
    if (st_use_electrode) then
       call af_add_cc_variable(tree, "lsf", ix=i_lsf)
    end if
 
    bc_method = "neumann_zero"
    call cfg_add_get(cfg, "species_boundary_condition", &
         bc_method, &
         "Boundary condition for the plasma species")
    select case (bc_method)
    case ("neumann_zero")
       bc_species => af_bc_neumann_zero
    case ("dirichlet_zero")
       bc_species => bc_species_dirichlet_zero
    case default
       print *, "Unknown boundary condition: ", trim(bc_method)
       print *, "Try neumann_zero or dirichlet_zero"
       error stop
    end select
 
    call cfg_add_get(cfg, "compute_power_density", compute_power_density, &
         "Whether to compute the deposited power density")
 
    if (compute_power_density) then
       call af_add_cc_variable(tree, "power_density", ix = i_power_density)
    end if
 
    call cfg_add_get(cfg, "use_end_streamer_length", st_use_end_streamer_length, &
         "Whether the length of the streamer is used to end the simulation")
    call cfg_add_get(cfg, "end_streamer_length", st_end_streamer_length, &
         "Streamer length at which the simulation will end.")
    call cfg_add_get(cfg, "initial_streamer_pos_steps_wait", &
         st_initial_streamer_pos_steps_wait, &
         "Number of simulation steps to wait before initializing "&
         "the starting position of the streamer")
 
    call cfg_add_get(cfg, "end_time", st_end_time, &
         "The desired endtime (s) of the simulation")
    call cfg_add_get(cfg, "box_size", st_box_size, &
         "The number of grid cells per coordinate in a box")
    call cfg_add_get(cfg, "coarse_grid_size", st_coarse_grid_size, &
         "The size of the coarse grid")
    call cfg_add_get(cfg, "domain_len", st_domain_len, &
         "The length of the domain (m)")
    call cfg_add_get(cfg, "domain_origin", st_domain_origin, &
         "The origin of the domain (m)")
    call cfg_add_get(cfg, "periodic", st_periodic, &
         "Whether the domain is periodic (per dimension)")
 
    call cfg_add_get(cfg, "plasma_region_enabled", st_plasma_region_enabled, &
         "Whether to limit plasma reactions to a certain region")
    call cfg_add_get(cfg, "plasma_region_rmin", st_plasma_region_rmin, &
         "Limit plasma reactions to coordinates between rmin and rmax")
    call cfg_add_get(cfg, "plasma_region_rmax", st_plasma_region_rmax, &
         "Limit plasma reactions to coordinates between rmin and rmax")
 
    if (all(st_coarse_grid_size == -1)) then
       ! Not set, automatically determine size
       st_coarse_grid_size = st_box_size * &
            nint(st_domain_len / minval(st_domain_len))
    end if
 
    tmp = maxval(st_domain_len/st_coarse_grid_size) / &
         minval(st_domain_len/st_coarse_grid_size)
    if (tmp > 1.001_dp) then
       print *, "!!! Warning: using non-square grid cells"
       write(*, "(A,F12.4)") " !!! Maximal aspect ratio:", tmp
    end if
 
    call cfg_add_get(cfg, "multigrid_num_vcycles", st_multigrid_num_vcycles, &
         "Number of V-cycles to perform per time step")
    call cfg_add_get(cfg, "multigrid_max_rel_residual", &
         st_multigrid_max_rel_residual, &
         "Stop multigrid when residual is smaller than this factor times max(|rhs|)")
 
    call cfg_add_get(cfg, "current_update_per_steps", &
         current_update_per_steps, &
         "Per how many iterations the electric current is computed")
 
    prolong_method = "limit"
    call cfg_add_get(cfg, "prolong_density", prolong_method, &
         "Density prolongation method (limit, linear, linear_cons, sparse)")
    select case (prolong_method)
    case ("limit")
       st_prolongation_method => af_prolong_limit
    case ("linear")
       st_prolongation_method => af_prolong_linear
    case ("linear_cons")
       st_prolongation_method => af_prolong_linear_cons
    case ("sparse")
       st_prolongation_method => af_prolong_sparse
    case ("zeroth")
       st_prolongation_method => af_prolong_zeroth
    case default
       error stop "Unknown prolong_density method"
    end select
 
    call cfg_add_get(cfg, "fixes%drt_max_field", st_drt_max_field, &
         "Enable flux limiting, but prevent field from exceeding this value")
    if (st_drt_max_field < 1e100_dp) then
       error stop "fixes%drt_max_field not yet implemented"
       st_drt_limit_flux = .true.
    end if
 
    call cfg_add_get(cfg, "fixes%source_factor", source_factor, &
         "Use source factor to prevent unphysical effects due to diffusion")
    call cfg_add_get(cfg, "fixes%write_source_factor", write_source_factor, &
         "Whether to write the source factor to the output")
    call cfg_add_get(cfg, "fixes%source_min_electrons_per_cell", &
         st_source_min_electrons_per_cell, &
         "Minimum number of electrons per cell to include source terms")
 
    select case (source_factor)
    case ("none")
       st_source_factor = source_factor_none
    case ("flux")
       st_source_factor = source_factor_flux
    case default
       print *, "Options fixes%source_factor: none, flux"
       error stop "Unknown fixes%source_factor"
    end select
 
    if (st_source_factor > 0 .and. write_source_factor) then
       call af_add_cc_variable(tree, "srcfac", ix=i_srcfac)
    end if
 
    call cfg_add_get(cfg, "rng_seed", rng_int4_seed, &
         "Seed for random numbers; if all zero, generate randomly")
 
    if (all(rng_int4_seed == 0)) then
       rng_int4_seed = get_random_seed()
       print *, "RNG seed: ", rng_int4_seed
 
       ! Store the updated seed in the configuration
       write(tmp_str, *) "rng_seed = ", rng_int4_seed
       call cfg_update_from_line(cfg, tmp_str)
    end if
 
    rng_int8_seed = transfer(rng_int4_seed, rng_int8_seed)
    call st_rng%set_seed(rng_int8_seed)
    n_threads = af_get_max_threads()
    call st_prng%init_parallel(n_threads, st_rng)
 
    ! Initialize global storage of reaction rates, +32 to avoid threads writing
    ! to nearby memory
    allocate(st_current_rates(n_reactions+32, n_threads))
    allocate(st_global_rates(n_reactions))
    allocate(st_current_jdote(1+32, n_threads))
    st_global_rates = 0
    st_global_jdote = 0
 
  end subroutine st_initialize
 
  !> Get a random seed based
  function get_random_seed() result(seed)
    use iso_fortran_env, only: int64
    integer  :: seed(4)
    integer  :: i
    real(dp) :: rr
    integer(int64) :: time
 
    ! Set a random seed (this does not work on all systems)
    call random_seed()
 
    ! Get some count of the time
    call system_clock(time)
 
    do i = 1, 4
       call random_number(rr)
       seed(i) = ieor(int(time), int(huge(1) * rr))
    end do
  end function get_random_seed
 
  !> Impose a Dirichlet zero boundary condition for plasma species in the last
  !> dimension, which is supposed to have the electrodes. We use Neumann
  !> conditions in the other dimensions. Note that this version avoids
  !> extrapolation (in contrast to the regular Dirichlet b.c.), which is more
  !> suitable for conserved species densities.
  subroutine bc_species_dirichlet_zero(box, nb, iv, coords, bc_val, bc_type)
    type(box_t), intent(in) :: box
    integer, intent(in)     :: nb
    integer, intent(in)     :: iv
    real(dp), intent(in)    :: coords(ndim, box%n_cell**(ndim-1))
    real(dp), intent(out)   :: bc_val(box%n_cell**(ndim-1))
    integer, intent(out)    :: bc_type
 
    if (af_neighb_dim(nb) == ndim) then
       bc_type = af_bc_dirichlet_copy
       bc_val  = 0.0_dp
    else
       bc_type = af_bc_neumann
       bc_val  = 0.0_dp
    end if
  end subroutine bc_species_dirichlet_zero
 
end module m_streamer