particles_to_grid_2d.f90

program particles_to_grid_2d
  use m_a2_all

  implicit none

  integer, parameter    :: box_size    = 8
  integer, parameter    :: i_phi       = 1
  integer, parameter    :: n_particles = 10*1000*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 :: coordinates(:, :), weights(:)

  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 particles_to_grid_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_adjust_refinement(tree, refine_routine, refine_info, 0)
     if (refine_info%n_add == 0) exit
  end do

  allocate(coordinates(2, n_particles))
  allocate(weights(n_particles))

  call random_number(coordinates(:, :))
  weights(:) = 1.0_dp

  do n = 1, n_particles
     coordinates(:, n) = (coordinates(:, n) - 0.5_dp) * domain_len
     ! coordinates(:, n) = (sqrt(coordinates(:, n)) - 0.5_dp) * domain_len
  end do

  call system_clock(t_start, count_rate)
  call a2_particles_to_grid(tree, i_phi, coordinates, weights, n_particles, 0)
  call system_clock(t_end, count_rate)
  call a2_write_silo(tree, "particles_to_grid_2d_0", 1, 0.0_dp, dir="output")
  call a2_tree_sum_cc(tree, i_phi, sum_density)
  print *, "zeroth order: ", (t_end-t_start) / real(count_rate, dp), " seconds"
  print *, "integrated density: ", sum_density

  call a2_tree_clear_cc(tree, i_phi)

  call system_clock(t_start, count_rate)
  call a2_particles_to_grid(tree, i_phi, coordinates, weights, n_particles, 1)
  call system_clock(t_end, count_rate)
  call a2_write_silo(tree, "particles_to_grid_2d_1", 1, 0.0_dp, dir="output")
  call a2_tree_sum_cc(tree, i_phi, sum_density)
  print *, "first order: ", (t_end-t_start) / real(count_rate, dp), " seconds"
  print *, "integrated density: ", sum_density

contains

  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 particles_to_grid_2d