solid_body_rotation.f90

#include "../src/cpp_macros.h"
 
program solid_body_rotation
  use m_af_all
  implicit none
 
  integer, parameter :: ncells = 8
  integer :: iq, i_flux
  integer, parameter :: coord_type = af_xyz
  real(dp) :: l_max(NDIM), l_min(NDIM)
  integer :: grid(NDIM)
  logical :: periodicBC(NDIM)
 
  type(af_t) :: tree
  !Add the stuff for AMR later
  real(dp) :: dt, time, end_time
  integer :: time_steps, t_iter
  character(len=100) :: fname
  print *, "Running solid body rotation 2D"
  print *, "Number of threads", af_get_max_threads()
 
  !Setting up the domain and ncells
  grid(:) = 8*ncells
  l_max(:) = 1.0_dp
  l_min(:) = -1.0_dp
  periodicbc(:) = .false.
 
 
  !Initialize the variables, boundary conditions and the tree
  call af_add_cc_variable(tree, "q", ix=iq)
  call af_add_fc_variable(tree, "flux", ix=i_flux)
  call af_set_cc_methods(tree, iq, af_bc_dirichlet_zero)
 
  call af_init(tree, ncells, l_max, grid, &
       periodic = periodicbc, r_min = l_min, &
       coord = coord_type)
 
  !Setting the boundary conditions
  call af_loop_box(tree, setinitconds)
 
  !Fill the ghost cells according to the initial conditions
  call af_gc_tree(tree, [iq])
 
  !Setting the time step data
  dt = 0.005_dp
  time = 0.0
  end_time = acos(-1.0_dp)
  time_steps = ceiling(end_time/dt)
 
  !Starting the time integration
  do t_iter = 1, time_steps
 
     if (mod(t_iter, 10) == 0) then
        write(fname, "(A,I0)") "output/solid_body_rotation_" // &
             dimname // "_", t_iter
        call af_write_silo(tree, trim(fname), t_iter, time)
     end if
 
 
     call af_loop_tree(tree, korenflux)
     !print *, 'Finished Looping over tree'
 
     call af_loop_box_arg(tree, updatesoln, [dt])
     !print *, 'finished time stepping'
 
     call af_gc_tree(tree, [iq])
     !print *, 'Updating ghost cells'
 
     time = time + dt
 
 
  end do
  call af_destroy(tree)
 
contains
 
  subroutine setinitconds( box )
    type(box_t), intent(inout) :: box
    integer :: IJK, nc
    real(dp) :: rr(NDIM)
 
    nc = box%n_cell
    do kji_do(0, nc+1)
       rr = af_r_cc(box, [ijk])
       box%cc(ijk, iq) = solution(rr)
    end do; close_do
  end subroutine setinitconds
 
  function solution( rr ) result(sol)
    real(dp), intent(in) :: rr(NDIM)
    real(dp) :: sol
    real(dp) :: peak(NDIM), xLim(NDIM), yLim(NDIM)
 
    peak(1) = 0.45_dp
    peak(2) = 0.0_dp
 
    xlim(1) = 0.1_dp
    xlim(2) = 0.6_dp
 
    ylim(1) = -0.25_dp
    ylim(2) = 0.25_dp
 
    if( xlim(1) < rr(1) .and. rr(1) < xlim(2) .and. ylim(1) < rr(2) &
         .and. rr(2) < ylim(2) ) then
       sol = 1.0_dp
    else if (norm2(rr + peak) < 0.35_dp) then
       sol = 1.0_dp - (norm2(rr + peak)/0.35_dp)
    else
       sol = 0.0_dp
    end if
 
  end function solution
 
  subroutine korenflux( tree, id )
    use m_af_flux_schemes
    type(af_t), intent(inout)  :: tree
    integer, intent(in) :: id
    integer :: nc, IJK
    real(dp), allocatable :: cc(DTIMES(:), :)
    real(dp), allocatable :: v(DTIMES(:),:)
    real(dp) :: rr(NDIM)
 
    nc = tree%boxes(id)%n_cell
    allocate(cc(dtimes(-1:nc+2), 1))
    allocate(v(dtimes(1:nc+1), ndim))
 
    !print *, "Inside Koren flux", id
    call af_gc2_box(tree, id, [iq], cc)
    !Computing the spatially varying velocities
    do kji_do(1,nc+1)
       rr = af_r_cc(tree%boxes(id), [ijk]) - 0.5_dp*tree%boxes(id)%dr
       v(ijk, 1) = 2.0_dp*rr(2)
       v(ijk, 2) = -2.0_dp*rr(1)
    end do; close_do
 
    call flux_koren_2d(cc(dtimes(:), iq), v, nc, 2)
    tree%boxes(id)%fc(:,:,:,i_flux) = v
 
  end subroutine korenflux
 
  subroutine updatesoln( box, dt )
    type(box_t), intent(inout) :: box
    real(dp), intent(in) :: dt(:)
    real(dp) :: inv_dr(NDIM)
    integer :: IJK, nc
 
    nc = box%n_cell
    inv_dr = 1/box%dr
 
    do j=1,nc
       do i=1,nc
          box%cc(i,j,iq) = box%cc(i,j,iq) - dt(1) * ( &
               inv_dr(1)* &
               (box%fc(i+1,j,1,iq) - box%fc(i,j,1,iq)) &
               + inv_dr(2)* &
               (box%fc(i,j+1,2,iq) - box%fc(i,j,2,iq)))
       end do
    end do
 
  end subroutine updatesoln
 
end program solid_body_rotation