WORKSPACE: Make workspace thread safe.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>

sources/amun.F90

@@ -171,7 +171,7 @@ program amun
     call print_message(loc, "Could not initialize module SYSTEM!")
   if (check_status(status /= 0)) go to 3000
 
-  call initialize_workspace(nwork, status)
+  call initialize_workspace(nwork, nthreads, status)
   if (status /= 0) &
     call print_message(loc, "Could not initialize module WORKSPACE!")
   if (check_status(status /= 0)) go to 2000

sources/evolution.F90

@@ -3526,6 +3526,9 @@ module evolution
 
     real(kind=8), dimension(NDIMS) :: dh, dhi
 
+    integer :: nt = 0
+!$  integer :: omp_get_thread_num
+
     real(kind=8), dimension(:,:,:,:,:)  , pointer, save :: f
     real(kind=8), dimension(:,:,:,:,:,:), pointer, save :: s
 
@@ -3533,6 +3536,7 @@ module evolution
 
 !-------------------------------------------------------------------------------
 !
+!$  nt = omp_get_thread_num
     if (first) then
       i = NDIMS * nf * nn**NDIMS
       j = 3 * i
@@ -3544,21 +3548,21 @@ module evolution
       end if
 
 #if NDIMS == 3
-      f(1:nf,1:nn,1:nn,1:nn,1:3)     => work(  1:i)
+      f(1:nf,1:nn,1:nn,1:nn,1:3)     => work(  1:i,nt)
-      s(1:nf,1:nn,1:nn,1:nn,1:2,1:3) => work(i+1:j)
+      s(1:nf,1:nn,1:nn,1:nn,1:2,1:3) => work(i+1:j,nt)
 #else /* NDIMS == 3 */
-      f(1:nf,1:nn,1:nn,1: 1,1:2)     => work(  1:i)
+      f(1:nf,1:nn,1:nn,1: 1,1:2)     => work(  1:i,nt)
-      s(1:nf,1:nn,1:nn,1: 1,1:2,1:2) => work(i+1:j)
+      s(1:nf,1:nn,1:nn,1: 1,1:2,1:2) => work(i+1:j,nt)
 #endif /* NDIMS == 3 */
 
       first = .false.
     end if
 
-    if (work_in_use) &
+    if (work_in_use(nt)) &
       call print_message(loc, "Workspace is being used right now! " //         &
                               "Corruptions can occur!")
 
-    work_in_use = .true.
+    work_in_use(nt) = .true.
 
     dh(1)  = adx(pdata%meta%level)
     dh(2)  = ady(pdata%meta%level)
@@ -3684,7 +3688,7 @@ module evolution
 #endif /* NDIMS == 3 */
     end if
 
-    work_in_use = .false.
+    work_in_use(nt) = .false.
 
     100 continue
 

sources/forcing.F90

@@ -1398,10 +1398,14 @@ module forcing
     real(kind=8), dimension(:,:,:,:), pointer, save :: acc
     real(kind=8), dimension(:,:,:)  , pointer, save :: den
 
+    integer :: nt = 0
+!$  integer :: omp_get_thread_num
+
     character(len=*), parameter :: loc = 'FORCING:inject_fmodes_block()'
 
 !-------------------------------------------------------------------------------
 !
+!$  nt = omp_get_thread_num
     if (first) then
       i = 3 * nn**NDIMS
       j = 4 * nn**NDIMS
@@ -1413,11 +1417,11 @@ module forcing
       end if
 
 #if NDIMS == 3
-      acc(1:3,1:nn,1:nn,1:nn) => work(  1:i)
+      acc(1:3,1:nn,1:nn,1:nn) => work(  1:i,nt)
-      den(1:nn,1:nn,1:nn)     => work(i+1:j)
+      den(1:nn,1:nn,1:nn)     => work(i+1:j,nt)
 #else /* NDIMS == 3 */
-      acc(1:2,1:nn,1:nn,1: 1) => work(1:i)
+      acc(1:2,1:nn,1:nn,1: 1) => work(  1:i,nt)
-      den(1:nn,1:nn,1: 1)     => work(i+1:j)
+      den(1:nn,1:nn,1: 1)     => work(i+1:j,nt)
 #endif /* NDIMS == 3 */
 
       first = .false.
@@ -1432,11 +1436,11 @@ module forcing
 #endif /* NDIMS == 3 */
     dvol = advol(pdata%meta%level)
 
-    if (work_in_use) &
+    if (work_in_use(nt)) &
       call print_message(loc, "Workspace is being used right now! " //         &
                               "Corruptions can occur!")
 
-    work_in_use = .true.
+    work_in_use(nt) = .true.
 
     acc(1:NDIMS,:,:,:) = 0.0d+00
 
@@ -1524,7 +1528,7 @@ module forcing
       pdata%u(ien,:,:,:) = pdata%u(ien,:,:,:) + den(:,:,:)
     end if
 
-    work_in_use = .false.
+    work_in_use(nt) = .false.
 
     100 continue
 

sources/mesh.F90

@@ -1901,12 +1901,16 @@ module mesh
     integer     , dimension(NDIMS) :: l, u
     real(kind=8), dimension(NDIMS) :: du
 
+    integer :: nt = 0
+!$  integer :: omp_get_thread_num
+
     real(kind=8), dimension(:,:,:), pointer, save :: tmp
 
     character(len=*), parameter :: loc = 'MESH::prolong_block()'
 
 !-------------------------------------------------------------------------------
 !
+!$  nt = omp_get_thread_num
     status = 0
 
     if (first) then
@@ -1918,16 +1922,16 @@ module mesh
         go to 100
       end if
 
-      tmp(1:pm(1),1:pm(2),1:pm(3)) => work(1:n)
+      tmp(1:pm(1),1:pm(2),1:pm(3)) => work(1:n,nt)
 
       first = .false.
     end if
 
-    if (work_in_use) &
+    if (work_in_use(nt)) &
       call print_message(loc, "Workspace is being used right now! " //         &
                               "Corruptions can occur!")
 
-    work_in_use = .true.
+    work_in_use(nt) = .true.
 
 #if NDIMS == 2
     k  = 1
@@ -2031,7 +2035,7 @@ module mesh
       end do ! nchildren
     end do ! n = 1, nv
 
-    work_in_use = .false.
+    work_in_use(nt) = .false.
 
     100 continue
 

sources/refinement.F90

@@ -508,10 +508,14 @@ module refinement
 
     real(kind=8), dimension(:,:,:,:), pointer, save :: w
 
+    integer :: nt = 0
+!$  integer :: omp_get_thread_num
+
     character(len=*), parameter :: loc = 'REFINEMENT:vorticity_magnitude()'
 
 !-------------------------------------------------------------------------------
 !
+!$  nt = omp_get_thread_num
     wmax = 0.0e+00
 
     if (first) then
@@ -524,9 +528,9 @@ module refinement
       end if
 
 #if NDIMS == 3
-      w(1:3,1:nn,1:nn,1:nn) => work(1:i)
+      w(1:3,1:nn,1:nn,1:nn) => work(1:i,nt)
 #else /* NDIMS == 3 */
-      w(1:3,1:nn,1:nn,1: 1) => work(1:i)
+      w(1:3,1:nn,1:nn,1: 1) => work(1:i,nt)
 #endif /* NDIMS == 3 */
 
       dh(:) = 1.0d+00
@@ -534,11 +538,11 @@ module refinement
       first = .false.
     end if
 
-    if (work_in_use) &
+    if (work_in_use(nt)) &
       call print_message(loc, &
                     "Workspace is being used right now! Corruptions can occur!")
 
-    work_in_use = .true.
+    work_in_use(nt) = .true.
 
     call curl(dh(:), pdata%q(ivx:ivz,:,:,:), w(:,:,:,:))
 
@@ -558,7 +562,7 @@ module refinement
     end do
 #endif /* NDIMS == 3 */
 
-    work_in_use = .false.
+    work_in_use(nt) = .false.
 
     wmax = sqrt(wmax)
 
@@ -606,10 +610,14 @@ module refinement
 
     real(kind=8), dimension(:,:,:,:), pointer, save :: w
 
+    integer :: nt = 0
+!$  integer :: omp_get_thread_num
+
     character(len=*), parameter :: loc = 'REFINEMENT:current_density_magnitude()'
 
 !-------------------------------------------------------------------------------
 !
+!$  nt = omp_get_thread_num
     jmax = 0.0e+00
 
     if (.not. magnetized) return
@@ -624,9 +632,9 @@ module refinement
       end if
 
 #if NDIMS == 3
-      w(1:3,1:nn,1:nn,1:nn) => work(1:i)
+      w(1:3,1:nn,1:nn,1:nn) => work(1:i,nt)
 #else /* NDIMS == 3 */
-      w(1:3,1:nn,1:nn,1: 1) => work(1:i)
+      w(1:3,1:nn,1:nn,1: 1) => work(1:i,nt)
 #endif /* NDIMS == 3 */
 
       dh(:) = 1.0d+00
@@ -634,11 +642,11 @@ module refinement
       first = .false.
     end if
 
-    if (work_in_use) &
+    if (work_in_use(nt)) &
       call print_message(loc, &
                     "Workspace is being used right now! Corruptions can occur!")
 
-    work_in_use = .true.
+    work_in_use(nt) = .true.
 
     call curl(dh(:), pdata%q(ibx:ibz,:,:,:), w(:,:,:,:))
 
@@ -658,7 +666,7 @@ module refinement
     end do
 #endif /* NDIMS == 3 */
 
-    work_in_use = .false.
+    work_in_use(nt) = .false.
 
     jmax = sqrt(jmax)
 

sources/sources.F90

@@ -296,6 +296,9 @@ module sources
     real(kind=8)  :: dvydx, dvydy, dvydz
     real(kind=8)  :: dvzdx, dvzdy, dvzdz
 
+    integer :: nt = 0
+!$  integer :: omp_get_thread_num
+
     real(kind=8), dimension(3)  :: ga, dh
     real(kind=8), dimension(nn) :: x, y
 #if NDIMS == 3
@@ -310,6 +313,7 @@ module sources
 
 !-------------------------------------------------------------------------------
 !
+!$  nt = omp_get_thread_num
     if (first) then
       i = nn**NDIMS
       j = 10 * i
@@ -321,11 +325,11 @@ module sources
       end if
 
 #if NDIMS == 3
-        db(1:nn,1:nn,1:nn)          => work(  1:i)
+        db(1:nn,1:nn,1:nn)          => work(  1:i,nt)
-        tmp(1:3,1:3,1:nn,1:nn,1:nn) => work(i+1:j)
+        tmp(1:3,1:3,1:nn,1:nn,1:nn) => work(i+1:j,nt)
 #else /* NDIMS == 3 */
-        db(1:nn,1:nn,1: 1)          => work(  1:i)
+        db(1:nn,1:nn,1: 1)          => work(  1:i,nt)
-        tmp(1:3,1:3,1:nn,1:nn,1: 1) => work(i+1:j)
+        tmp(1:3,1:3,1:nn,1:nn,1: 1) => work(i+1:j,nt)
 #endif /* NDIMS == 3 */
 
       first = .false.
@@ -399,11 +403,11 @@ module sources
 !
     if (viscosity > 0.0d+00) then
 
-      if (work_in_use) &
+      if (work_in_use(nt)) &
         call print_message(loc, "Workspace is being used right now! " //       &
                                 "Corruptions can occur!")
 
-      work_in_use = .true.
+      work_in_use(nt) = .true.
 
 ! prepare coordinate increments
 !
@@ -531,7 +535,7 @@ module sources
 
       end if ! ien > 0
 
-      work_in_use = .false.
+      work_in_use(nt) = .false.
 
     end if ! viscosity is not zero
 
@@ -539,11 +543,11 @@ module sources
 !
     if (magnetized) then
 
-      if (work_in_use) &
+      if (work_in_use(nt)) &
         call print_message(loc, "Workspace is being used right now! " //       &
                                 "Corruptions can occur!")
 
-      work_in_use = .true.
+      work_in_use(nt) = .true.
 
 ! prepare coordinate increments
 !
@@ -718,7 +722,7 @@ module sources
 
       end if ! resistivity is not zero
 
-      work_in_use = .false.
+      work_in_use(nt) = .false.
 
     end if ! magnetized
 

sources/statistics.F90

@@ -421,10 +421,14 @@ module statistics
     integer(kind=4), dimension(nprocs) :: cdist
 #endif /* MPI */
 
+    integer :: nt = 0
+!$  integer :: omp_get_thread_num
+
     character(len=*), parameter :: loc = 'INTEGRALS:store_statistics()'
 
 !-------------------------------------------------------------------------------
 !
+!$  nt = omp_get_thread_num
 ! process and store the mesh statistics only on the master node
 !
     if (master) then
@@ -490,16 +494,16 @@ module statistics
       n = ni**NDIMS
       l = 1
       u = n
-      vel(1:ni,1:ni,1:nk) => work(l:u)
+      vel(1:ni,1:ni,1:nk) => work(l:u,nt)
       l = l + n
       u = u + n
-      mag(1:ni,1:ni,1:nk) => work(l:u)
+      mag(1:ni,1:ni,1:nk) => work(l:u,nt)
       l = l + n
       u = u + n
-      sqd(1:ni,1:ni,1:nk) => work(l:u)
+      sqd(1:ni,1:ni,1:nk) => work(l:u,nt)
       l = l + n
       u = u + n
-      tmp(1:ni,1:ni,1:nk) => work(l:u)
+      tmp(1:ni,1:ni,1:nk) => work(l:u,nt)
 
       first = .false.
     end if
@@ -526,11 +530,11 @@ module statistics
       mxarr(7) = 0.0d+00
     end if
 
-    if (work_in_use) &
+    if (work_in_use(nt)) &
       call print_message(loc, "Workspace is being used right now! " //         &
                               "Corruptions can occur!")
 
-    work_in_use = .true.
+    work_in_use(nt) = .true.
 
 ! associate the pointer with the first block on the data block list
 !
@@ -699,7 +703,7 @@ module statistics
 
     end do ! data blocks
 
-    work_in_use = .false.
+    work_in_use(nt) = .false.
 
 #ifdef MPI
 ! sum the integral array from all processes

sources/workspace.F90

@@ -33,13 +33,17 @@ module workspace
 !
   integer, save :: nwork = 0
 
+! the last thread number
+!
+  integer, save :: nt = 0
+
 ! the flag indicating that the workspace is in use
 !
-  logical, save :: work_in_use = .false.
+  logical, dimension(:), allocatable, save :: work_in_use
 
 ! the common workspace to use for local arrays
 !
-  real(kind=8), dimension(:), allocatable, target :: work
+  real(kind=8), dimension(:,:), allocatable, target :: work
 
   private
 
@@ -64,18 +68,19 @@ module workspace
 !
 !   Arguments:
 !
-!     ninit  - the initial workspace size;
+!     ninit    - the initial workspace size;
-!     status - the call status (0 for success, otherwise failure);
+!     nthreads - the number of threads;
+!     status   - the call status (0 for success, otherwise failure);
 !
 !===============================================================================
 !
-  subroutine initialize_workspace(ninit, status)
+  subroutine initialize_workspace(ninit, nthreads, status)
 
     use helpers, only : print_message
 
     implicit none
 
-    integer, intent(in)  :: ninit
+    integer, intent(in)  :: ninit, nthreads
     integer, intent(out) :: status
 
     character(len=*), parameter :: loc = 'WORKSPACE::initialize_workspace()'
@@ -91,12 +96,15 @@ module workspace
     end if
 
     nwork = ninit
+    nt    = nthreads - 1
 
-    allocate(work(nwork), stat=status)
+    allocate(work_in_use(0:nt), work(nwork,0:nt), stat=status)
 
     if (status /= 0) &
       call print_message(loc, "Could not allocate the common workspace!")
 
+    work_in_use = .false.
+
 !-------------------------------------------------------------------------------
 !
   end subroutine initialize_workspace
@@ -127,7 +135,7 @@ module workspace
     status = 0
 
     if (allocated(work)) then
-      deallocate(work, stat=status)
+      deallocate(work_in_use, work, stat=status)
 
       if (status /= 0) &
         call print_message(loc, "Could not deallocate the common workspace!")
@@ -164,7 +172,7 @@ module workspace
 !
     status = 0
 
-    if (work_in_use) then
+    if (any(work_in_use)) then
       call print_message(loc, "Could not resize the workspace. " //            &
                               "It is being used right now!")
 
@@ -173,7 +181,7 @@ module workspace
       if (nsize > nwork) then
         deallocate(work, stat=status)
         if (status == 0) then
-          allocate(work(nsize), stat=status)
+          allocate(work(nsize,0:nt), stat=status)
           if (status /= 0) then
             call print_message(loc, "Could not allocate a new workspace!")
             status = 1