m_refine.f90

Go to the documentation of this file.

!> Module with grid refinement settings and routines
module m_refine
#include "../afivo/src/cpp_macros.h"
  use m_af_all
 
  implicit none
  private
 
  ! The refinement buffer width in cells (around flagged cells)
  integer, public, protected :: refine_buffer_width = 4
 
  ! The number of steps after which the mesh is updated
  integer, public, protected :: refine_per_steps = 2
 
  ! The grid spacing will always be larger than this value
  real(dp), public, protected :: refine_min_dx = 1.0e-7_dp
 
  ! The grid spacing will always be smaller than this value
  real(dp), public, protected :: refine_max_dx = 1.0e-3_dp
 
  ! Refine if alpha*dx is larger than this value
  real(dp), protected :: refine_adx = 1.0_dp
 
  ! For refinement, use alpha(f * E)/f, where f is this factor
  real(dp), protected :: refine_adx_fac = 1.0_dp
 
  ! Refine if the curvature in phi is larger than this value
  real(dp), protected :: refine_cphi = 1e99_dp
 
  ! Allow derefinement if the curvature in phi is smaller than this value
  real(dp), protected :: derefine_cphi = 1e99_dp
 
  ! Only derefine if grid spacing if smaller than this value
  real(dp), protected :: derefine_dx = 1e-4_dp
 
  ! Refine around initial conditions up to this time
  real(dp), protected :: refine_init_time = 10e-9_dp
 
  ! Refine until dx is smaller than this factor times the seed width
  real(dp), protected :: refine_init_fac = 0.25_dp
 
  ! Ensure grid spacing around electrode is less than this value
  real(dp), public :: refine_electrode_dx = 1e99_dp
 
  ! Auxiliary variable used for electrode derefinement
  ! it will be the same value as above
  real(dp), public :: current_electrode_dx = 1e99_dp
 
  ! Multiply the above value by this factor to derefine 
  ! during interpulse
  real(dp), public, protected :: electrode_derefine_factor = 1.0
 
  ! Start refining electrode some time before the start of a new pulse
  real(dp), public, protected :: refine_prepulse_time = 1.0e-9_dp
 
  ! Minimum electron density for adding grid refinement
  real(dp), protected :: refine_min_dens = -1.0e99_dp
 
  ! Refine a region up to this grid spacing
  real(dp), protected, allocatable :: refine_regions_dr(:)
 
  ! Refine regions up to this simulation time
  real(dp), protected, allocatable :: refine_regions_tstop(:)
 
  ! Minimum coordinate of the refinement regions
  real(dp), protected, allocatable :: refine_regions_rmin(:,:)
 
  ! Maximum coordinate of the refinement regions
  real(dp), protected, allocatable :: refine_regions_rmax(:,:)
 
  ! Limit refinement in a region to this grid spacing
  real(dp), protected, allocatable :: refine_limits_dr(:)
 
  ! Minimum coordinate of the refinement limits
  real(dp), protected, allocatable :: refine_limits_rmin(:,:)
 
  ! Maximum coordinate of the refinement limits
  real(dp), protected, allocatable :: refine_limits_rmax(:,:)
 
  ! Use effective alpha (minus attachment) for refinement
  logical, protected :: refine_use_alpha_effective
 
  procedure(af_subr_ref), pointer :: refine_routine => null()
 
  ! Public methods
  public :: refine_initialize
  public :: refine_routine
 
contains
 
  !> Initialize the grid refinement options
  subroutine refine_initialize(cfg)
    use m_config
    use m_gas
    use m_user_methods
    type(cfg_t), intent(inout) :: cfg
    integer                    :: n
    real(dp)                   :: vec(ndim)
    real(dp), allocatable      :: dbuffer(:)
 
    !< [basic_refinement_settings]
    call cfg_add_get(cfg, "refine_buffer_width", refine_buffer_width, &
         "The refinement buffer width in cells (around flagged cells)")
    call cfg_add_get(cfg, "refine_per_steps", refine_per_steps, &
         "The number of steps after which the mesh is updated")
    call cfg_add_get(cfg, "refine_min_dx", refine_min_dx, &
         "The grid spacing will always be larger than this value (m)")
    call cfg_add_get(cfg, "refine_max_dx", refine_max_dx, &
         "The grid spacing will always be smaller than this value (m)")
    call cfg_add_get(cfg, "refine_adx", refine_adx, &
         "Refine if alpha*dx is larger than this value")
    call cfg_add_get(cfg, "derefine_dx", derefine_dx, &
         "Only derefine if grid spacing if smaller than this value")
    call cfg_add_get(cfg, "refine_init_time", refine_init_time, &
         "Refine around initial conditions up to this time")
    call cfg_add_get(cfg, "refine_init_fac", refine_init_fac, &
         "Refine until dx is smaller than this factor times the seed width")
    call cfg_add_get(cfg, "refine_electrode_dx", refine_electrode_dx, &
         "Ensure grid spacing around electrode is less than this value (m)")
    !< [basic_refinement_settings]
 
    if (refine_min_dx > refine_max_dx) &
         error stop "Cannot have refine_min_dx < refine_max_dx"
    current_electrode_dx = refine_electrode_dx
 
 
    call cfg_add_get(cfg, "refine_adx_fac", refine_adx_fac, &
         "For refinement, use alpha(f * E)/f, where f is this factor")
    call cfg_add_get(cfg, "refine_cphi", refine_cphi, &
         "Refine if the curvature in phi is larger than this value")
    call cfg_add_get(cfg, "derefine_cphi", derefine_cphi, &
         "Allow derefinement if the curvature in phi is smaller than this value")
 
    call cfg_add_get(cfg, "electrode_derefine_factor", &
         electrode_derefine_factor, &
         "Multiplication factor to derefine electrode during interpulse")
    call cfg_add_get(cfg, "refine_prepulse_time", refine_prepulse_time, &
         "Start refining electrode some time before the start of a new pulse")
    call cfg_add_get(cfg, "refine_min_dens", refine_min_dens, &
         "Minimum electron density for adding grid refinement")
    call cfg_add_get(cfg, "refine_use_alpha_effective", refine_use_alpha_effective, &
         "Use effective alpha (minus attachment) for refinement")
 
    call cfg_add(cfg, "refine_regions_dr", [1.0e99_dp], &
         "Refine regions up to this grid spacing (m)", .true.)
    call cfg_add(cfg, "refine_regions_tstop", [1.0e99_dp], &
         "Refine regions up to this simulation time", .true.)
    vec = 0.0_dp
    call cfg_add(cfg, "refine_regions_rmin", vec, &
         "Minimum coordinate of the refinement regions", .true.)
    call cfg_add(cfg, "refine_regions_rmax", vec, &
         "Maximum coordinate of the refinement regions", .true.)
 
    call cfg_get_size(cfg, "refine_regions_dr", n)
    allocate(refine_regions_dr(n))
    allocate(refine_regions_tstop(n))
    allocate(refine_regions_rmin(ndim, n))
    allocate(refine_regions_rmax(ndim, n))
    allocate(dbuffer(ndim * n))
 
    call cfg_get(cfg, "refine_regions_dr", refine_regions_dr)
    call cfg_get(cfg, "refine_regions_tstop", refine_regions_tstop)
    call cfg_get(cfg, "refine_regions_rmin", dbuffer)
    refine_regions_rmin = reshape(dbuffer, [ndim, n])
    call cfg_get(cfg, "refine_regions_rmax", dbuffer)
    refine_regions_rmax = reshape(dbuffer, [ndim, n])
 
    call cfg_add(cfg, "refine_limits_dr", [1.0e99_dp], &
         "Refine regions at most up to this grid spacing", .true.)
    vec = 0.0_dp
    call cfg_add(cfg, "refine_limits_rmin", vec, &
         "Minimum coordinate of the refinement limits", .true.)
    call cfg_add(cfg, "refine_limits_rmax", vec, &
         "Maximum coordinate of the refinement limits", .true.)
 
    call cfg_get_size(cfg, "refine_limits_dr", n)
    allocate(refine_limits_dr(n))
    allocate(refine_limits_rmin(ndim, n))
    allocate(refine_limits_rmax(ndim, n))
    deallocate(dbuffer)
    allocate(dbuffer(ndim * n))
 
    call cfg_get(cfg, "refine_limits_dr", refine_limits_dr)
    call cfg_get(cfg, "refine_limits_rmin", dbuffer)
    refine_limits_rmin = reshape(dbuffer, [ndim, n])
    call cfg_get(cfg, "refine_limits_rmax", dbuffer)
    refine_limits_rmax = reshape(dbuffer, [ndim, n])
 
    if (associated(user_refine)) then
       refine_routine => user_refine
    else
       refine_routine => default_refinement
    end if
 
  end subroutine refine_initialize
 
  !> Set the cell refinement flags for box
  subroutine default_refinement(box, cell_flags)
    use m_streamer
    use m_geometry
    use m_init_cond
    use m_gas
    use m_lookup_table
    use m_transport_data
    type(box_t), intent(in) :: box
    !> Refinement flags for the cells of the box
    integer, intent(out)    :: &
         cell_flags(dtimes(box%n_cell))
    integer                 :: ijk, n, nc
    real(dp)                :: min_dx, max_dx, gas_dens
    real(dp)                :: alpha, adx, fld, elec_dens
    real(dp)                :: dist, rmin(ndim), rmax(ndim)
 
    nc = box%n_cell
    min_dx = minval(box%dr)
    max_dx = maxval(box%dr)
 
    do kji_do(1,nc)
       if (gas_constant_density) then
          gas_dens = gas_number_density
       else
          gas_dens = box%cc(ijk, i_gas_dens)
       end if
       fld   = box%cc(ijk, i_electric_fld) * si_to_townsend / gas_dens
 
       if (refine_use_alpha_effective) then
          alpha = (lt_get_col(td_tbl, td_alpha, refine_adx_fac * fld) - &
               lt_get_col(td_tbl, td_eta, refine_adx_fac * fld)) * &
               gas_dens / refine_adx_fac
          alpha = max(0.0_dp, alpha)
       else
          alpha = lt_get_col(td_tbl, td_alpha, refine_adx_fac * fld) * &
               gas_dens / refine_adx_fac
       end if
 
       adx   = max_dx * alpha
       elec_dens = box%cc(ijk, i_electron)
 
       if (adx > refine_adx .and. elec_dens > refine_min_dens) then
          cell_flags(ijk) = af_do_ref
       else if (adx < 0.125_dp * refine_adx .and. max_dx < derefine_dx) then
          cell_flags(ijk) = af_rm_ref
       else
          cell_flags(ijk) = af_keep_ref
       end if
 
       ! Refine around the initial conditions
       if (global_time < refine_init_time) then
          do n = 1, init_conds%n_cond
             dist = gm_dist_line(af_r_cc(box, [ijk]), &
                  init_conds%seed_r0(:, n), &
                  init_conds%seed_r1(:, n), ndim)
             if (dist - init_conds%seed_width(n) < 2 * max_dx &
                  .and. max_dx > refine_init_fac * &
                  init_conds%seed_width(n)) then
                cell_flags(ijk) = af_do_ref
             end if
          end do
       end if
 
       ! Refine around electrode
       if (iand(box%tag, mg_lsf_box) > 0 .and. &
            max_dx > current_electrode_dx) then
          cell_flags(ijk) = af_do_ref
       end if
    end do; close_do
 
    ! Check fixed refinements
    rmin = box%r_min
    rmax = box%r_min + box%dr * box%n_cell
 
    do n = 1, size(refine_regions_dr)
       if (global_time <= refine_regions_tstop(n) .and. &
            max_dx > refine_regions_dr(n) .and. all(&
            rmax >= refine_regions_rmin(:, n) .and. &
            rmin <= refine_regions_rmax(:, n))) then
          ! Mark just the center cell to prevent refining neighbors
          cell_flags(dtimes(nc/2)) = af_do_ref
       end if
    end do
 
    do n = 1, size(refine_limits_dr)
       if (max_dx < 2 * refine_limits_dr(n) .and. all(&
            rmin >= refine_limits_rmin(:, n) .and. &
            rmax <= refine_limits_rmax(:, n))) then
          ! Mark just the center cell to prevent refining neighbors
          where (cell_flags == af_do_ref) cell_flags = af_keep_ref
       end if
    end do
 
    ! Make sure we don't have or get a too fine or too coarse grid
    if (max_dx > refine_max_dx) then
       cell_flags = af_do_ref
    else if (min_dx < 2 * refine_min_dx) then
       where(cell_flags == af_do_ref) cell_flags = af_keep_ref
    end if
 
  end subroutine default_refinement
 
end module m_refine