diff --git a/sources/amun.F90 b/sources/amun.F90 index 725714e..fa8506d 100644 --- a/sources/amun.F90 +++ b/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 diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 0359800..76fbbf3 100644 --- a/sources/evolution.F90 +++ b/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) - s(1:nf,1:nn,1:nn,1:nn,1:2,1:3) => work(i+1:j) + 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,nt) #else /* NDIMS == 3 */ - f(1:nf,1:nn,1:nn,1: 1,1:2) => work( 1:i) - s(1:nf,1:nn,1:nn,1: 1,1:2,1:2) => work(i+1:j) + 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,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 diff --git a/sources/forcing.F90 b/sources/forcing.F90 index 1978ff8..e73abd0 100644 --- a/sources/forcing.F90 +++ b/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) - den(1:nn,1:nn,1:nn) => work(i+1:j) + acc(1:3,1:nn,1:nn,1:nn) => work( 1:i,nt) + 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) - den(1:nn,1:nn,1: 1) => work(i+1:j) + acc(1:2,1:nn,1:nn,1: 1) => work( 1:i,nt) + 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 diff --git a/sources/mesh.F90 b/sources/mesh.F90 index 32a00f7..e75e48b 100644 --- a/sources/mesh.F90 +++ b/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 diff --git a/sources/refinement.F90 b/sources/refinement.F90 index 44aaad3..984efcc 100644 --- a/sources/refinement.F90 +++ b/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) diff --git a/sources/sources.F90 b/sources/sources.F90 index a846941..67d03e1 100644 --- a/sources/sources.F90 +++ b/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) - tmp(1:3,1:3,1:nn,1:nn,1:nn) => work(i+1:j) + 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,nt) #else /* NDIMS == 3 */ - db(1:nn,1:nn,1: 1) => work( 1:i) - tmp(1:3,1:3,1:nn,1:nn,1: 1) => work(i+1:j) + 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,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 diff --git a/sources/statistics.F90 b/sources/statistics.F90 index c802bef..d0e8040 100644 --- a/sources/statistics.F90 +++ b/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 diff --git a/sources/workspace.F90 b/sources/workspace.F90 index 515fd15..46f7a8a 100644 --- a/sources/workspace.F90 +++ b/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; -! status - the call status (0 for success, otherwise failure); +! ninit - the initial workspace size; +! 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