blocking

README

@@ -89,6 +89,15 @@ list of selected runtime parameters and default values.
   domain decomposition is in y and z-directions; when it is 3, the
   domain decomposition is in 3D.
 
+* tb_blocksize_x = -1
+* tb_blocksize_y = 16
+* tb_blocksize_z = 16
+
+  SMC uses a blocking strategy.  These parameters determines blocking
+  size in x, y and z-directions.  The value should be either greater
+  than or equal to 4, or less than 0.  If it is less than 0, no
+  blocking is used in that direction.
+
 * verbose = 0
 
   This determines verbosity.

inputs_SMC

@@ -5,6 +5,11 @@
  n_cellz = 128
  max_grid_size = 128
 
+ ! block size inside a box
+ tb_blocksize_x = -1   ! no blocking in x-direction
+ tb_blocksize_y = 16
+ tb_blocksize_z = 16
+
  verbose = 2
 
  stop_time = 3.d-3

src/advance.f90

@@ -173,7 +173,7 @@ subroutine dUdt (U, Uprime, Q, mu, xi, lam, Ddiag, dx, courno, istep)
     double precision, intent(inout), optional :: courno
 
     integer :: lo(3), hi(3)
-    integer :: n
+    integer :: n, iblock
 
     logical :: update_courno
     double precision :: courno_proc
@@ -240,7 +240,7 @@ subroutine dUdt (U, Uprime, Q, mu, xi, lam, Ddiag, dx, courno, istep)
     !
     ! Hyperbolic and Transport terms
     !
-    !$omp parallel private(n,lo,hi,up,ulo,uhi,upp,uplo,uphi) &
+    !$omp parallel private(n,iblock,lo,hi,up,ulo,uhi,upp,uplo,uphi) &
     !$omp private(qp,qlo,qhi,mup,xip,lamp,Ddp)
     do n=1,nfabs(Q)
 
@@ -261,14 +261,17 @@ subroutine dUdt (U, Uprime, Q, mu, xi, lam, Ddiag, dx, courno, istep)
        uplo = lbound(upp)
        uphi = ubound(upp)
 
-       lo = tb_get_valid_lo(n)
-       hi = tb_get_valid_hi(n)
+       do iblock = 1, tb_get_nblocks(n)
+          lo = tb_get_block_lo(iblock,n)
+          hi = tb_get_block_hi(iblock,n)
 
-       call narrow_diffterm_3d(lo,hi,dx,qp,qlo(1:3),qhi(1:3),upp,uplo(1:3),uphi(1:3), &
-            mup,xip,lamp,Ddp)
+          call narrow_diffterm_3d(lo,hi,dx,qp,qlo(1:3),qhi(1:3),upp,uplo(1:3),uphi(1:3), &
+               mup,xip,lamp,Ddp)
 
-       call hypterm_3d(lo,hi,dx,up,ulo(1:3),uhi(1:3),qp,qlo(1:3),qhi(1:3),&
-            upp,uplo(1:3),uphi(1:3))
+          call hypterm_3d(lo,hi,dx,up,ulo(1:3),uhi(1:3),qp,qlo(1:3),qhi(1:3),&
+               upp,uplo(1:3),uphi(1:3))
+
+       end do
 
     end do
     !$omp end parallel

src/probin.f90

@@ -31,6 +31,9 @@ module probin_module
   logical, save, public :: tb_collapse_boxes = .false.
   integer, save, public :: tb_idim_more = 2
   integer, save, public :: tb_idim_less = 1
+  integer, save, public :: tb_blocksize_x = -1
+  integer, save, public :: tb_blocksize_y = 16
+  integer, save, public :: tb_blocksize_z = 16
 
   ! These will be allocated and defined below
   logical,    allocatable, save, public :: pmask(:)
@@ -72,6 +75,9 @@ module runtime_init_module
   namelist /probin/ tb_collapse_boxes
   namelist /probin/ tb_idim_more
   namelist /probin/ tb_idim_less
+  namelist /probin/ tb_blocksize_x
+  namelist /probin/ tb_blocksize_y
+  namelist /probin/ tb_blocksize_z
 
   private
 

src/threadbox.f90

@@ -16,19 +16,39 @@ module threadbox_module
   integer, allocatable, save :: tb_lo(:,:), tb_hi(:,:) 
   integer, allocatable, save :: tb_glo(:,:), tb_ghi(:,:) 
 
-  !$omp threadprivate(worktodo,tb_lo,tb_hi,tb_glo,tb_ghi)
+  type blocks  ! in_a_threadbox
+     integer :: nblocks
+     integer, pointer :: lo(:,:) => Null()
+     integer, pointer :: hi(:,:) => Null()
+  end type blocks
+
+  type(blocks), allocatable, save :: allblocks(:)
+
+  !$omp threadprivate(worktodo,tb_lo,tb_hi,tb_glo,tb_ghi,allblocks)
 
   private
 
   public ::build_threadbox, destroy_threadbox,  &
        tb_get_valid_lo, tb_get_valid_hi, tb_get_grown_lo, tb_get_grown_hi, &
+       tb_get_block_lo, tb_get_block_hi, tb_get_nblocks, &
        tb_multifab_setval, tb_worktodo
 
 contains
 
   subroutine destroy_threadbox()
     integer :: ibox
     !$omp parallel private(ibox)
+    if (allocated(allblocks)) then
+       do ibox=1,nb
+          if (associated(allblocks(ibox)%lo)) then
+             deallocate(allblocks(ibox)%lo)
+          end if
+          if (associated(allblocks(ibox)%hi)) then
+             deallocate(allblocks(ibox)%hi)
+          end if
+       end do
+       deallocate(allblocks)
+    end if
     if (allocated(tb_lo)) deallocate(tb_lo)
     if (allocated(tb_hi)) deallocate(tb_hi)
     if (allocated(tb_glo)) deallocate(tb_glo)
@@ -38,7 +58,8 @@ subroutine destroy_threadbox()
   end subroutine destroy_threadbox
 
   subroutine build_threadbox(la, ng_in)
-    use probin_module, only : tb_split_dim, tb_collapse_boxes, tb_idim_more, tb_idim_less
+    use probin_module, only : tb_split_dim, tb_collapse_boxes, tb_idim_more, tb_idim_less, &
+         tb_blocksize_x, tb_blocksize_y, tb_blocksize_z
     implicit none
     type(layout), intent(in) :: la
     integer, intent(in) :: ng_in
@@ -74,6 +95,7 @@ subroutine build_threadbox(la, ng_in)
     allocate(tb_glo(3,nb))
     allocate(tb_ghi(3,nb))
     allocate(worktodo(nb))
+    allocate(allblocks(nb))
 
     !$omp single
     call setup_boxgroups(tb_collapse_boxes) 
@@ -92,6 +114,8 @@ subroutine build_threadbox(la, ng_in)
 
     call init_threadbox(la)
 
+    call init_allblocks(tb_blocksize_x,tb_blocksize_y,tb_blocksize_z)
+
     !$omp end parallel
 
   end subroutine build_threadbox
@@ -370,6 +394,82 @@ subroutine split_domain(totalsize, n, isize, start)
 
   end subroutine split_domain
 
+
+  subroutine init_allblocks(blksizex,blksizey,blksizez)
+    implicit none
+    integer,intent(in) :: blksizex,blksizey,blksizez
+    integer :: idim, ibox, nbk(3), tbsize(3)
+    integer :: i,j,k,ibk,nblk
+    integer :: bksize(3), zero_lo(3), zero_hi(3)
+    integer, allocatable :: xsize(:), ysize(:), zsize(:)
+    integer, allocatable :: xstart(:), ystart(:), zstart(:)
+
+    bksize(1) = blksizex
+    bksize(2) = blksizey
+    bksize(3) = blksizez
+
+    do ibox=1,nb
+
+       if (.not. worktodo(ibox)) then
+
+          allblocks(ibox)%nblocks = 0
+
+       else
+
+          tbsize = tb_hi(:,ibox) - tb_lo(:,ibox) + 1
+
+          do idim=1,3
+             if (bksize(idim) <= 0) then
+                nbk(idim) = 1
+             else
+                nbk(idim) = max(int(tbsize(idim)/bksize(idim)),1)
+             end if
+          end do
+
+          nblk = nbk(1)*nbk(2)*nbk(3)
+          allblocks(ibox)%nblocks = nblk
+          allocate(allblocks(ibox)%lo(3,nblk))
+          allocate(allblocks(ibox)%hi(3,nblk))
+
+          allocate(xsize (nbk(1)))
+          allocate(ysize (nbk(2)))
+          allocate(zsize (nbk(3)))
+          allocate(xstart(nbk(1)))
+          allocate(ystart(nbk(2)))
+          allocate(zstart(nbk(3)))
+
+          call split_domain(tbsize(1), nbk(1), xsize, xstart)
+          call split_domain(tbsize(2), nbk(2), ysize, ystart)
+          call split_domain(tbsize(3), nbk(3), zsize, zstart)
+
+          ibk = 1
+          do k = 1, nbk(3)
+             do j = 1, nbk(2)
+                do i = 1, nbk(1)
+
+                   zero_lo(1) = xstart(i)
+                   zero_lo(2) = ystart(j)
+                   zero_lo(3) = zstart(k)
+                   zero_hi(1) = xstart(i) + xsize(i) - 1
+                   zero_hi(2) = ystart(j) + ysize(j) - 1
+                   zero_hi(3) = zstart(k) + zsize(k) - 1
+
+                   allblocks(ibox)%lo(:,ibk) = zero_lo + tb_lo(:,ibox)
+                   allblocks(ibox)%hi(:,ibk) = zero_hi + tb_lo(:,ibox)
+
+                   ibk = ibk + 1
+                end do
+             end do
+          end do
+
+          deallocate(xsize,ysize,zsize,xstart,ystart,zstart)
+
+       end if
+    end do
+
+  end subroutine init_allblocks
+
+
   function tb_get_valid_lo(ilocal) result (lo)
     implicit none
     integer, intent(in) :: ilocal
@@ -400,6 +500,28 @@ function tb_get_grown_hi(ilocal) result (hi)
   end function tb_get_grown_hi
 
 
+  function tb_get_nblocks(ilocal) result (nblk)
+    implicit none
+    integer, intent(in) :: ilocal
+    integer :: nblk
+    nblk = allblocks(ilocal)%nblocks
+  end function tb_get_nblocks
+
+  function tb_get_block_lo(iblock, ilocal) result(lo)
+    implicit none
+    integer, intent(in) :: iblock, ilocal
+    integer, dimension(3) :: lo
+    lo = allblocks(ilocal)%lo(:,iblock)
+  end function tb_get_block_lo
+
+  function tb_get_block_hi(iblock, ilocal) result(hi)
+    implicit none
+    integer, intent(in) :: iblock, ilocal
+    integer, dimension(3) :: hi
+    hi = allblocks(ilocal)%hi(:,iblock)
+  end function tb_get_block_hi
+
+
   subroutine tb_multifab_setval(mf, val, all)
     type(multifab), intent(inout) :: mf
     real(dp_t), intent(in) :: val