photon_absorption_2d.f90

program photon_absorption_2d
  use m_a2_all

  implicit none

  integer, parameter    :: box_size   = 8
  integer, parameter    :: i_phi      = 1
  integer, parameter    :: n_photons  = 1000
  real(dp), parameter   :: domain_len = 1.0_dp
  real(dp), parameter   :: r_min(2)  = -0.5_dp * domain_len
  real(dp), parameter   :: dr         = domain_len / box_size
  real(dp), allocatable :: coord_t0(:, :), coord_t1(:, :)

  type(a2_t)      :: tree
  type(ref_info_t) :: refine_info
  integer          :: n, id, ix_list(2, 1)
  integer          :: count_rate, t_start, t_end
  real(dp)         :: sum_density

  print *, "Running photon_absorption_2d"
  print *, "Number of threads", af_get_max_threads()

  ! Initialize tree
  call a2_init(tree, & ! Tree to initialize
       box_size, &     ! A box contains box_size**DIM cells
       n_var_cell=1, & ! Number of cell-centered variables
       n_var_face=0, & ! Number of face-centered variables
       dr=dr, &        ! Distance between cells on base level
       cc_names=["phi"], & ! Variable names
       r_min=r_min)

  ! Set up geometry
  id             = 1
  ix_list(:, id) = 1            ! Set index of box

  ! Create the base mesh, using the box indices and their neighbor information
  call a2_set_base(tree, 1, ix_list)

  call a2_set_cc_methods(tree, i_phi, a2_bc_neumann_zero, a2_gc_interp)

  do
     call a2_loop_box(tree, set_opaque_cells)
     call a2_adjust_refinement(tree, refine_routine, refine_info, 0)
     if (refine_info%n_add == 0) exit
  end do

  ! Get photon start and end position in domain
  allocate(coord_t0(2, n_photons))
  allocate(coord_t1(2, n_photons))

  call random_number(coord_t0(:, :))
  call random_number(coord_t1(:, :))

  coord_t0 = (coord_t0 - 0.5_dp) * domain_len
  coord_t1 = (coord_t1 - 0.5_dp) * domain_len

  call system_clock(t_start, count_rate)

  call system_clock(t_end, count_rate)
  call a2_write_silo(tree, "photon_absorption_2d_0", 1, 0.0_dp, dir="output")

contains

  subroutine set_opaque_cells(box)
    type(box2_t), intent(inout) :: box
    integer                      :: i, j, nc
    real(dp)                     :: rr(2)

    nc = box%n_cell

    do j = 0, nc+1; do i = 0, nc+1
       rr = a2_r_cc(box, [i, j])

       ! Change phi to a 'high' value in part of the domain
       if (norm2(rr) < 0.25_dp) then
          box%cc(i, j, i_phi) = 2.0_dp
       else
          box%cc(i, j, i_phi) = 1.0_dp
       end if
    end do; end do

  end subroutine set_opaque_cells

  subroutine refine_routine(box, cell_flags)
    type(box2_t), intent(in) :: box
    integer, intent(out)      :: cell_flags(box%n_cell, box%n_cell)

    if (all(box%r_min < 0.0_dp) .and. box%lvl <= 5) then
       cell_flags(:, :) = af_do_ref
    else
       cell_flags(:, :) = af_keep_ref
    end if
  end subroutine refine_routine

end program photon_absorption_2d