use threadbox

GNUmakefile

@@ -13,7 +13,7 @@ COMP := GNU
 
 ### 
 ifdef MPI
-  # Set USE_MPI_WRAPPERS := t to use mpif90 or mpiifort, 
+  # Set USE_MPI_WRAPPERS := t to use mpif90 or mpiifot, 
   # otherwise you need to specify mpi_include_dir, mpi_lib_dir, and mpi_libraries.
   USE_MPI_WRAPPERS := t
   ifndef USE_MPI_WRAPPERS

README

@@ -70,6 +70,15 @@ list of selected runtime parameters and default values.
   and the number of MPI tasks.  Ideally, you would like to have one
   box for each MPI task to minimize the communication cost.
 
+* tb_split_dim = 2
+
+  This determines how domain decomposition in some subroutines is done
+  for threads.  When OpenMP is used, each box is "virtually" divided
+  into "nthreads" boxes and each OpenMP thread works on one
+  thread-box.  This parameter can be either 2 or 3. When it is 2, the
+  domain decomposition is in y and z-directions; when it is 3, the
+  domain decomposition is in 3D.
+
 * verbose = 0
 
   This determines verbosity.

src/advance.f90

@@ -3,6 +3,7 @@ module advance_module
   use bl_error_module
   use kernels_module
   use multifab_module
+  use threadbox_module
   use time_module
   use transport_properties
   use variables_module
@@ -27,7 +28,7 @@ subroutine advance(U, dt, courno, dx, istep)
     double precision,  intent(in   ) :: dx(3)
     integer, intent(in) :: istep
 
-    call multifab_setval(Unew, 0.d0, .true.)
+    call tb_multifab_setval(Unew, 0.d0, .true.)
 
     call dUdt(U, Uprime, Q, mu, xi, lam, Ddiag, dx, courno, istep)
     call set_dt(dt, courno, istep)
@@ -130,16 +131,18 @@ subroutine update_rk3 (a,U1,b,U2,c,Uprime)
 
     nc = ncomp(U1)
 
+    !$omp parallel private(i,j,k,m,n,lo,hi,u1p,u2p,upp)
     do n=1,nfabs(U1)
 
+       if (.not.tb_worktodo(n)) cycle
+
        u1p => dataptr(U1,    n)
        u2p => dataptr(U2,    n)
        upp => dataptr(Uprime,n)
 
-       lo = lwb(get_box(Uprime,n))
-       hi = upb(get_box(Uprime,n))
+       lo = tb_get_valid_lo(n)
+       hi = tb_get_valid_hi(n)
 
-       !$omp parallel do private(i,j,k,m) collapse(3)
        do m = 1, nc
           do k = lo(3),hi(3)
              do j = lo(2),hi(2)
@@ -149,8 +152,8 @@ subroutine update_rk3 (a,U1,b,U2,c,Uprime)
              end do
           end do
        end do
-       !$omp end parallel do
     end do
+    !$omp end parallel
 
   end subroutine update_rk3
 
@@ -191,7 +194,7 @@ subroutine dUdt (U, Uprime, Q, mu, xi, lam, Ddiag, dx, courno, istep)
     call multifab_fill_boundary(U)
     wt_fillboundary = wt_fillboundary + (parallel_wtime()-wt1)
 
-    call multifab_setval(Uprime, 0.d0)
+    call tb_multifab_setval(Uprime, 0.d0)
 
     !
     ! Calculate primitive variables based on U
@@ -237,8 +240,12 @@ 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 private(qp,qlo,qhi,mup,xip,lamp,Ddp)
     do n=1,nfabs(Q)
 
+       if (.not.tb_worktodo(n)) cycle
+
        up => dataptr(U,n)
        upp=> dataptr(Uprime,n)
        qp => dataptr(Q,n)
@@ -254,8 +261,8 @@ subroutine dUdt (U, Uprime, Q, mu, xi, lam, Ddiag, dx, courno, istep)
        uplo = lbound(upp)
        uphi = ubound(upp)
 
-       lo = lwb(get_box(Q,n))
-       hi = upb(get_box(Q,n))
+       lo = tb_get_valid_lo(n)
+       hi = tb_get_valid_hi(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)
@@ -264,6 +271,7 @@ subroutine dUdt (U, Uprime, Q, mu, xi, lam, Ddiag, dx, courno, istep)
             upp,uplo(1:3),uphi(1:3))
 
     end do
+    !$omp end parallel
     wt2 = parallel_wtime()
     wt_hypdiff = wt_hypdiff + (wt2-wt1)
 
@@ -281,18 +289,23 @@ subroutine compute_courno(Q, dx, courno)
     integer :: n, lo(3), hi(3), qlo(4), qhi(4)
     double precision, pointer :: qp(:,:,:,:)
 
+    !$omp parallel private(n, lo, hi, qlo, qhi, qp) &
+    !$omp reduction(max:courno)
     do n=1,nfabs(Q)
 
+       if (.not.tb_worktodo(n)) cycle
+
        qp => dataptr(Q,n)
        qlo = lbound(qp)
        qhi = ubound(qp)
 
-       lo = lwb(get_box(Q,n))
-       hi = upb(get_box(Q,n))
+       lo = tb_get_valid_lo(n)
+       hi = tb_get_valid_hi(n)
 
        call comp_courno_3d(lo,hi,dx,qp,qlo(1:3),qhi(1:3),courno)
 
     end do
+    !$omp end parallel
   end subroutine compute_courno
 
 end module advance_module

src/initialize.f90

@@ -2,6 +2,7 @@ module initialize_module
 
   use multifab_module
   use variables_module
+  use threadbox_module
 
   implicit none
 
@@ -52,6 +53,8 @@ subroutine initialize_from_scratch(la,dt,courno,dx,U)
 
     ng = stencil_ng
 
+    call build_threadbox(la,ng)
+
     call multifab_build(U,la,ncons,ng)
 
     call init_data(U,dx,prob_lo,prob_hi)