random_refinement.f90

This example shows how to create an AMR tree, perform random refinement, write output files, and how to fill ghost cells.

#include "../src/cpp_macros.h"
!> \example random_refinement.f90
!>
!> This example shows how to create an AMR tree, perform random refinement,
!> write output files, and how to fill ghost cells.
program random_refinement
  use m_af_all
 
  implicit none
 
  type(af_t)         :: tree
  integer            :: grid_size(NDIM)
  real(dp)           :: domain_size(NDIM)
  logical            :: periodic(NDIM) = .true.
  integer            :: iter, boxes_used
  integer, parameter :: coord_type     = af_xyz ! af_xyz or af_cyl
 
  integer, parameter :: box_size   = 8
  integer, parameter :: i_phi      = 1
  type(ref_info_t)   :: ref_info
  real(dp)           :: sum_phi_t0, sum_phi
  character(len=100) :: fname
  integer            :: count_rate,t_start,t_end
 
  write(*,'(A,I0,A)') 'program random_refinement_', ndim, "d"
 
  print *, "Number of threads", af_get_max_threads()
 
  call af_add_cc_variable(tree, "phi")
 
  call af_set_cc_methods(tree, 1, af_bc_dirichlet_zero, &
       prolong=af_prolong_linear)
 
  ! Initialize tree
  grid_size(1) = 2 * box_size
  grid_size(2:) = box_size
  domain_size(1) = 2 * acos(-1.0_dp)
  domain_size(2:) = acos(-1.0_dp)
 
  if (command_argument_count() == 1) then
     ! Load tree from file
     call get_command_argument(1, fname)
     call af_read_tree(tree, fname, read_other_data=read_string)
  else
     call af_init(tree, & ! Tree to initialize
          box_size, &     ! A box contains box_size**DIM cells
          domain_size, &
          grid_size, &
          periodic=periodic, &
          coord=coord_type)
 
     call af_print_info(tree)
 
     ! Set variables on base by using the helper functions af_loop_box(tree, sub)
     ! and af_loop_boxes(tree, sub). These functions call the subroutine sub for
     ! each box in the tree, with a slightly different syntax.
     call af_loop_box(tree, set_init_cond)
 
     ! Fill ghost cells for phi. The third argument is a subroutine that fills
     ! ghost cells near refinement boundaries, and the fourth argument fill ghost
     ! cells near physical boundaries.
     call af_gc_tree(tree, [i_phi])
  end if
 
  call af_tree_sum_cc(tree, i_phi, sum_phi_t0)
 
  call system_clock(t_start, count_rate)
  boxes_used = 1
  do iter = 1, 50
     ! This writes a VTK output file containing the cell-centered values of the
     ! leaves of the tree (the boxes not covered by refinement).
     ! Variables are the names given as the third argument.
     write(fname, "(A,I0)") "output/random_refinement_" // dimname // "_", iter
     call af_write_silo(tree, trim(fname), n_cycle=iter)
     call af_write_tree(tree, fname, write_other_data=write_string)
 
     ! This updates the refinement of the tree, by at most one level per call.
     ! The second argument is a subroutine that is called for each box that can
     ! be refined or derefined, and it should set refinement flags. Information
     ! about the changes in refinement are returned in the third argument.
     !
     ! Within each (de)refinement step the subroutine af_tidy_up is called. This
     ! means that every now and then whether too much holes appear in the tree
     ! box list. Therefore reordering of tree is not necessary at the end of the
     ! iteration process
     call af_adjust_refinement(tree, ref_routine, ref_info, ref_buffer=0)
 
     call af_tree_sum_cc(tree, i_phi, sum_phi)
     write(*, "(A,E10.2)") " conservation error: ", sum_phi - sum_phi_t0
     boxes_used = boxes_used + ref_info%n_add - ref_info%n_rm
     write(*,'(4(3x,A,1x,i6))') "# new     boxes", ref_info%n_add, &
          "# removed boxes", ref_info%n_rm,  &
          "# boxes used   ", boxes_used,     &
          " highest level ", tree%highest_lvl
     print *, "boxes / in use / not used: ", tree%highest_id, &
          count(tree%boxes(1:tree%highest_id)%in_use), &
          tree%n_removed_ids
  end do
  call system_clock(t_end, count_rate)
 
  write(*, '(A,i3,1x,A,f8.2,1x,A,/)') &
       ' Wall-clock time after ',iter, &
       ' iterations: ', (t_end-t_start) / real(count_rate, dp), &
       ' seconds'
 
  call af_print_info(tree)
 
  ! This call is not really necessary here, but cleaning up the data in a tree
  ! is important if your program continues with other tasks.
  call af_destroy(tree)
 
contains
 
  ! Return the refinement flag for boxes(id)
  subroutine ref_routine(box, cell_flags)
    type(box_t), intent(in) :: box ! A list of all boxes in the tree
    integer, intent(out)     :: cell_flags(DTIMES(box%n_cell))
    real(dp)                 :: rr
 
    ! Draw a [0, 1) random number
    call random_number(rr)
 
    if (rr < 0.5_dp**ndim .and. box%lvl < 5) then
       cell_flags = af_do_ref   ! Add refinement
    else
       cell_flags = af_rm_ref   ! Ask to remove this box, which will not always
       ! happen (see documentation)
    end if
  end subroutine ref_routine
 
  ! This routine sets the initial conditions for each box
  subroutine set_init_cond(box)
    type(box_t), intent(inout) :: box
    integer                     :: IJK, nc
    real(dp)                    :: rr(NDIM)
 
    nc = box%n_cell
    do kji_do(1,nc)
       ! Get the coordinate of the cell center at i,j
       rr = af_r_cc(box, [ijk])
 
       ! Set the values at each cell according to some function
#if NDIM > 1
          box%cc(ijk, i_phi) = sin(0.5_dp * rr(1))**2 * cos(rr(2))**2
#else
          box%cc(ijk, i_phi) = sin(0.5_dp * rr(1))**2
#endif
    end do; close_do
  end subroutine set_init_cond
 
  subroutine write_string(my_unit)
    integer, intent(in) :: my_unit
    character(len=10)   :: str
 
    str = "Hello"
    write(my_unit) str
  end subroutine write_string
 
  subroutine read_string(my_unit)
    integer, intent(in) :: my_unit
    character(len=10)   :: str
    read(my_unit) str
 
    if (str /= "Hello") error stop "error reading other data"
  end subroutine read_string
 
end program random_refinement