Revert b296377

Change back to using threadprivate for xyzcache
ADACS-Australia · Mar 11, 2022 · 82a3284 · 82a3284
1 parent 25e4434
commit 82a3284
Show file tree

Hide file tree

Showing 5 changed files with 30 additions and 85 deletions.
diff --git a/build/Makefile b/build/Makefile
@@ -1781,7 +1781,7 @@ endif
 SOURCES= physcon.f90 ${CONFIG} ${SRCKERNEL} io.F90 units.f90 boundary.f90 \
          mpi_utils.F90 dtype_kdtree.F90 utils_omp.F90 utils_cpuinfo.f90 \
          utils_allocate.f90 icosahedron.f90 \
-         utils_mathfunc.f90 part.F90 omp_cache.f90 ${DOMAIN} utils_timing.f90 mpi_balance.F90 \
+         utils_mathfunc.f90 part.F90 ${DOMAIN} utils_timing.f90 mpi_balance.F90 \
          setup_params.f90 timestep.f90 utils_dumpfiles.f90 utils_indtimesteps.F90 utils_infiles.f90 \
          utils_sort.f90 utils_supertimestep.F90 utils_tables.f90 utils_gravwave.f90 \
          utils_sphNG.f90 utils_vectors.f90 utils_datafiles.f90 datafiles.f90 \
@@ -1861,9 +1861,9 @@ edit: checksetup
 # these are the sources for anything which uses the readwrite_dumps module
 #
 SRCDUMP= physcon.f90 ${CONFIG} ${SRCKERNEL} io.F90 units.f90 boundary.f90 mpi_utils.F90 \
-         utils_timing.f90 utils_infiles.f90 dtype_kdtree.f90 utils_allocate.f90 part.F90 \
-         omp_cache.f90 ${DOMAIN} mpi_dens.F90 mpi_force.F90 \
-         mpi_balance.F90 mpi_memory.F90 mpi_derivs.F90 kdtree.F90 linklist_kdtree.F90 \
+         utils_timing.f90 utils_infiles.f90 dtype_kdtree.f90 utils_allocate.f90 part.F90 ${DOMAIN} \
+         mpi_dens.F90 mpi_force.F90 \
+         mpi_balance.F90 mpi_stack.F90 mpi_derivs.F90 kdtree.F90 linklist_kdtree.F90 \
          utils_dumpfiles.f90 utils_vectors.f90 utils_mathfunc.f90 \
          utils_datafiles.f90 utils_filenames.f90 utils_tables.f90 datafiles.f90 gitinfo.f90 \
          centreofmass.f90 \
@@ -2301,7 +2301,7 @@ cleanphantomsinks:
 # these are the sources for the multirun utility
 #
 SRCMULT = physcon.f90 ${CONFIG} ${SRCKERNEL} io.F90 mpi_utils.F90 utils_timing.f90 dtype_kdtree.F90 ${SRCFASTMATH} \
-          units.f90 boundary.f90 utils_allocate.f90 part.F90 omp_cache.f90 timestep.f90 setup_params.f90 \
+          units.f90 boundary.f90 utils_allocate.f90 part.F90 timestep.f90 setup_params.f90 \
           utils_filenames.f90 utils_mathfunc.f90 utils_vectors.f90 utils_omp.F90 utils_datafiles.f90 datafiles.f90 utils_tables.f90 utils_infiles.f90 \
           ${SRCEOS} viscosity.f90 options.f90 damping.f90 utils_gravwave.f90 \
           utils_dumpfiles.f90 utils_summary.f90 centreofmass.f90 \

diff --git a/src/main/dens.F90 b/src/main/dens.F90
@@ -99,7 +99,7 @@ module densityforce
 
  !--kernel related parameters
  !real, parameter    :: cnormk = 1./pi, wab0 = 1., gradh0 = -3.*wab0, radkern2 = 4F.0
-
+ integer, parameter :: isizecellcache = 50000
  integer, parameter :: isizeneighcache = 0
  integer, parameter :: maxdensits = 50
 
@@ -137,8 +137,6 @@ subroutine densityiterate(icall,npart,nactive,xyzh,vxyzu,divcurlv,divcurlB,Bevol
  use part,      only:ngradh
  use viscosity, only:irealvisc
  use io_summary,only:summary_variable,iosumhup,iosumhdn
- use omp_cache, only:ompcache,maxcellcache
-
  integer,      intent(in)    :: icall,npart,nactive
  real,         intent(inout) :: xyzh(:,:)
  real,         intent(in)    :: vxyzu(:,:),fxyzu(:,:),fext(:,:)
@@ -152,6 +150,9 @@ subroutine densityiterate(icall,npart,nactive,xyzh,vxyzu,divcurlv,divcurlB,Bevol
  real,         intent(inout) :: radprop(:,:)
  real(kind=4), intent(out)   :: dvdx(:,:)
 
+ real,   save :: xyzcache(isizecellcache,3)
+!$omp threadprivate(xyzcache)
+
  integer :: i,icell
  integer :: nneigh,np,npcell
  integer :: nwarnup,nwarndown,nwarnroundoff
@@ -173,16 +174,14 @@ subroutine densityiterate(icall,npart,nactive,xyzh,vxyzu,divcurlv,divcurlB,Bevol
  integer                   :: n_remote_its,nlocal
  real                      :: ntotal
  logical                   :: iterations_finished,do_export
- integer                   :: ithread,OMP_GET_THREAD_NUM
-
 
  if (mpi) then
     call init_cell_exchange(xrecvbuf,irequestrecv)
     call reset_stacks
  endif
 
  if (iverbose >= 3 .and. id==master) &
-    write(iprint,*) ' cell cache =',maxcellcache,' neigh cache = ',isizeneighcache,' icall = ',icall
+    write(iprint,*) ' cell cache =',isizecellcache,' neigh cache = ',isizeneighcache,' icall = ',icall
 
  if (icall==0 .or. icall==1) then
     call reset_neighbour_stats(nneightry,nneighact,maxneightry,maxneighact,ncalc,nrelink)
@@ -284,11 +283,7 @@ subroutine densityiterate(icall,npart,nactive,xyzh,vxyzu,divcurlv,divcurlB,Bevol
 !$omp reduction(+:nrelink) &
 !$omp reduction(+:stressmax) &
 !$omp reduction(max:rhomax) &
-!$omp shared(ompcache) &
-!$omp private(i,ithread)
-
- ithread = OMP_GET_THREAD_NUM() + 1
-
+!$omp private(i)
 
 !$omp do schedule(runtime)
  over_cells: do icell=1,int(ncells)
@@ -300,7 +295,7 @@ subroutine densityiterate(icall,npart,nactive,xyzh,vxyzu,divcurlv,divcurlB,Bevol
     !
     !--get the neighbour list and fill the cell cache
     !
-    call get_neighbour_list(icell,listneigh,nneigh,xyzh,ompcache(:,1:3,ithread),maxcellcache,getj=.false., &
+    call get_neighbour_list(icell,listneigh,nneigh,xyzh,xyzcache,isizecellcache,getj=.false., &
                            remote_export=remote_export)
     do_export = any(remote_export)
     cell%icell                   = icell
@@ -323,7 +318,7 @@ subroutine densityiterate(icall,npart,nactive,xyzh,vxyzu,divcurlv,divcurlB,Bevol
        !$omp end critical (send_and_recv_remote)
     endif
 
-    call compute_cell(cell,listneigh,nneigh,getdv,getdB,Bevol,xyzh,vxyzu,fxyzu,fext,ompcache(:,1:3,ithread),rad)
+    call compute_cell(cell,listneigh,nneigh,getdv,getdB,Bevol,xyzh,vxyzu,fxyzu,fext,xyzcache,rad)
 
     if (do_export) then
        ! write directly to stack
@@ -337,7 +332,7 @@ subroutine densityiterate(icall,npart,nactive,xyzh,vxyzu,divcurlv,divcurlB,Bevol
           if (.not. converged) then
              if (redo_neighbours) then
                 call set_hmaxcell(cell%icell,cell%hmax)
-                call get_neighbour_list(-1,listneigh,nneigh,xyzh,ompcache(:,1:3,ithread),maxcellcache,getj=.false., &
+                call get_neighbour_list(-1,listneigh,nneigh,xyzh,xyzcache,isizecellcache,getj=.false., &
                                       cell_xpos=cell%xpos,cell_xsizei=cell%xsizei,cell_rcuti=cell%rcuti, &
                                       remote_export=remote_export)
 
@@ -352,7 +347,7 @@ subroutine densityiterate(icall,npart,nactive,xyzh,vxyzu,divcurlv,divcurlB,Bevol
                 nrelink = nrelink + 1
              endif
 
-             call compute_cell(cell,listneigh,nneigh,getdv,getdB,Bevol,xyzh,vxyzu,fxyzu,fext,ompcache(:,1:3,ithread),rad)
+             call compute_cell(cell,listneigh,nneigh,getdv,getdB,Bevol,xyzh,vxyzu,fxyzu,fext,xyzcache,rad)
 
              if (do_export) then
                 stack_waiting%cells(cell%waiting_index) = cell
@@ -405,10 +400,10 @@ subroutine densityiterate(icall,npart,nactive,xyzh,vxyzu,divcurlv,divcurlB,Bevol
           cell = stack_remote%cells(i)
 
           ! icell is unused (-1 here)
-          call get_neighbour_list(-1,listneigh,nneigh,xyzh,ompcache(:,1:3,ithread),maxcellcache,getj=.false., &
+          call get_neighbour_list(-1,listneigh,nneigh,xyzh,xyzcache,isizecellcache,getj=.false., &
                                   cell_xpos=cell%xpos,cell_xsizei=cell%xsizei,cell_rcuti=cell%rcuti)
 
-          call compute_cell(cell,listneigh,nneigh,getdv,getdB,Bevol,xyzh,vxyzu,fxyzu,fext,ompcache(:,1:3,ithread),rad)
+          call compute_cell(cell,listneigh,nneigh,getdv,getdB,Bevol,xyzh,vxyzu,fxyzu,fext,xyzcache,rad)
 
           remote_export = .false.
           remote_export(cell%owner+1) = .true. ! use remote_export array to send back to the owner
@@ -457,7 +452,7 @@ subroutine densityiterate(icall,npart,nactive,xyzh,vxyzu,divcurlv,divcurlB,Bevol
           !$omp end critical (send_and_recv_remote)
           if (.not. converged) then
              call set_hmaxcell(cell%icell,cell%hmax)
-             call get_neighbour_list(-1,listneigh,nneigh,xyzh,ompcache(:,1:3,ithread),maxcellcache,getj=.false., &
+             call get_neighbour_list(-1,listneigh,nneigh,xyzh,xyzcache,isizecellcache,getj=.false., &
                                     cell_xpos=cell%xpos,cell_xsizei=cell%xsizei,cell_rcuti=cell%rcuti, &
                                     remote_export=remote_export)
 
@@ -466,7 +461,7 @@ subroutine densityiterate(icall,npart,nactive,xyzh,vxyzu,divcurlv,divcurlB,Bevol
              call reserve_stack(stack_redo,cell%waiting_index)
              call send_cell(cell,remote_export,irequestsend,xsendbuf) ! send the cell to remote
              !$omp end critical (send_and_recv_remote)
-             call compute_cell(cell,listneigh,nneigh,getdv,getdB,Bevol,xyzh,vxyzu,fxyzu,fext,ompcache(:,1:3,ithread),rad)
+             call compute_cell(cell,listneigh,nneigh,getdv,getdB,Bevol,xyzh,vxyzu,fxyzu,fext,xyzcache,rad)
 
              stack_redo%cells(cell%waiting_index) = cell
           else
@@ -556,7 +551,6 @@ pure subroutine get_density_sums(i,xpartveci,hi,hi1,hi21,iamtypei,iamgasi,iamdus
  use kernel,   only:get_kernel,get_kernel_grav1
  use part,     only:iphase,iamgas,iamdust,iamtype,maxphase,ibasetype,igas,idust,rhoh,massoftype,iradxi
  use dim,      only:ndivcurlv,gravity,maxp,nalpha,use_dust,do_radiation
- use omp_cache,only:maxcellcache
  integer,      intent(in)    :: i
  real,         intent(in)    :: xpartveci(:)
  real(kind=8), intent(in)    :: hi,hi1,hi21
@@ -621,7 +615,7 @@ pure subroutine get_density_sums(i,xpartveci,hi,hi1,hi21,iamtypei,iamgasi,iamdus
     if (ifilledneighcache .and. n <= isizeneighcache) then
        rij2 = dxcache(1,n)
     else
-       if (ifilledcellcache .and. n <= maxcellcache) then
+       if (ifilledcellcache .and. n <= isizecellcache) then
           ! positions from cache are already mod boundary
           dx = xpartveci(ixi) - xyzcache(n,1)
           dy = xpartveci(iyi) - xyzcache(n,2)
@@ -1155,7 +1149,6 @@ pure subroutine compute_cell(cell,listneigh,nneigh,getdv,getdB,Bevol,xyzh,vxyzu,
  use viscosity,   only:irealvisc
  use io,          only:id
  use dim,         only:mpi
- use omp_cache,   only:maxcellcache
 
  type(celldens),  intent(inout)  :: cell
 
@@ -1165,7 +1158,7 @@ pure subroutine compute_cell(cell,listneigh,nneigh,getdv,getdB,Bevol,xyzh,vxyzu,
  logical,         intent(in)     :: getdB
  real,            intent(in)     :: Bevol(:,:)
  real,            intent(in)     :: xyzh(:,:),vxyzu(:,:),fxyzu(:,:),fext(:,:)
- real,            intent(in)     :: xyzcache(maxcellcache,3)
+ real,            intent(in)     :: xyzcache(isizecellcache,3)
  real,            intent(in)     :: rad(:,:)
 
  real                            :: dxcache(7,isizeneighcache)

diff --git a/src/main/force.F90 b/src/main/force.F90
@@ -56,6 +56,8 @@ module forces
  character(len=80), parameter, public :: &  ! module version
     modid="$Id$"
 
+ integer, parameter :: maxcellcache = 50000
+
  public :: force, reconstruct_dv ! latter to avoid compiler warning
 
  !--indexing for xpartveci array
@@ -222,7 +224,6 @@ subroutine force(icall,npart,xyzh,vxyzu,fxyzu,divcurlv,divcurlB,Bevol,dBevol,&
  use mpimemory,    only:stack_remote  => force_stack_1
  use mpimemory,    only:stack_waiting => force_stack_2
  use io_summary,   only:iosumdtr
- use omp_cache,    only:ompcache,maxcellcache
  integer,      intent(in)    :: icall,npart
  real,         intent(in)    :: xyzh(:,:)
  real,         intent(inout) :: vxyzu(:,:)
@@ -242,6 +243,8 @@ subroutine force(icall,npart,xyzh,vxyzu,fxyzu,divcurlv,divcurlB,Bevol,dBevol,&
  real,         intent(inout) :: radprop(:,:)
  real,         intent(in)    :: dens(:), metrics(:,:,:,:)
 
+ real, save :: xyzcache(maxcellcache,4)
+!$omp threadprivate(xyzcache)
  integer :: i,icell,nneigh
  integer :: nstokes,nsuper,ndrag,ndustres
  real    :: dtmini,dtohm,dthall,dtambi,dtvisc
@@ -275,7 +278,6 @@ subroutine force(icall,npart,xyzh,vxyzu,fxyzu,divcurlv,divcurlB,Bevol,dBevol,&
  logical                   :: remote_export(nprocs), do_export
  type(cellforce)           :: cell,xrecvbuf(nprocs),xsendbuf
  integer                   :: irequestsend(nprocs),irequestrecv(nprocs)
- integer                   :: ithread,OMP_GET_THREAD_NUM
 
 #ifdef IND_TIMESTEPS
  nbinmaxnew      = 0
@@ -429,11 +431,7 @@ subroutine force(icall,npart,xyzh,vxyzu,fxyzu,divcurlv,divcurlB,Bevol,dBevol,&
 !$omp shared(dustfrac) &
 !$omp shared(ddustevol) &
 !$omp shared(deltav) &
-!$omp shared(ibin_wake,ibinnow_m1) &
-!$omp shared(ompcache) &
-!$omp private(ithread)
-
- ithread = OMP_GET_THREAD_NUM() + 1
+!$omp shared(ibin_wake,ibinnow_m1)
 
 !$omp do schedule(runtime)
  over_cells: do icell=1,int(ncells)
@@ -454,7 +452,7 @@ subroutine force(icall,npart,xyzh,vxyzu,fxyzu,divcurlv,divcurlB,Bevol,dBevol,&
     !--get the neighbour list and fill the cell cache
     !
 
-    call get_neighbour_list(icell,listneigh,nneigh,xyzh,ompcache(:,:,ithread),maxcellcache,getj=.true., &
+    call get_neighbour_list(icell,listneigh,nneigh,xyzh,xyzcache,maxcellcache,getj=.true., &
 #ifdef GRAVITY
                            f=cell%fgrav, &
 #endif
@@ -477,7 +475,7 @@ subroutine force(icall,npart,xyzh,vxyzu,fxyzu,divcurlv,divcurlB,Bevol,dBevol,&
 
     call compute_cell(cell,listneigh,nneigh,Bevol,xyzh,vxyzu,fxyzu, &
                       iphase,divcurlv,divcurlB,alphaind,eta_nimhd,eos_vars, &
-                      dustfrac,dustprop,gradh,ibinnow_m1,ibin_wake,stressmax,ompcache(:,:,ithread),&
+                      dustfrac,dustprop,gradh,ibinnow_m1,ibin_wake,stressmax,xyzcache,&
                       rad,radprop,dens,metrics)
 
     if (do_export) then
@@ -518,15 +516,15 @@ subroutine force(icall,npart,xyzh,vxyzu,fxyzu,divcurlv,divcurlB,Bevol,dBevol,&
     over_remote: do i = 1,stack_remote%n
        cell = stack_remote%cells(i)
 
-       call get_neighbour_list(-1,listneigh,nneigh,xyzh,ompcache(:,:,ithread),maxcellcache,getj=.true., &
+       call get_neighbour_list(-1,listneigh,nneigh,xyzh,xyzcache,maxcellcache,getj=.true., &
 #ifdef GRAVITY
                          f=cell%fgrav, &
 #endif
                          cell_xpos=cell%xpos,cell_xsizei=cell%xsizei,cell_rcuti=cell%rcuti)
 
        call compute_cell(cell,listneigh,nneigh,Bevol,xyzh,vxyzu,fxyzu, &
                          iphase,divcurlv,divcurlB,alphaind,eta_nimhd,eos_vars, &
-                         dustfrac,dustprop,gradh,ibinnow_m1,ibin_wake,stressmax,ompcache(:,:,ithread),&
+                         dustfrac,dustprop,gradh,ibinnow_m1,ibin_wake,stressmax,xyzcache,&
                          rad,radprop,dens,metrics)
 
        remote_export = .false.
@@ -865,7 +863,6 @@ subroutine compute_forces(i,iamgasi,iamdusti,xpartveci,hi,hi1,hi21,hi41,gradhi,g
 #endif
  use utils_gr,    only:get_bigv
  use radiation_utils, only:get_rad_R
- use omp_cache,   only:maxcellcache
  integer,         intent(in)    :: i
  logical,         intent(in)    :: iamgasi,iamdusti
  real,            intent(in)    :: xpartveci(:)

diff --git a/src/main/memory.F90 b/src/main/memory.F90
@@ -30,7 +30,6 @@ subroutine allocate_memory(n, part_only)
  use allocutils, only:nbytes_allocated,bytes2human
  use part,       only:allocate_part
  use linklist,   only:allocate_linklist,ifirstincell
- use omp_cache,  only:allocate_cache
  use mpimemory,  only:allocate_mpi_memory
  use mpibalance, only:allocate_balance_arrays
  use mpiderivs,  only:allocate_comms_arrays
@@ -93,8 +92,6 @@ subroutine allocate_memory(n, part_only)
     endif
  endif
 
- call allocate_cache
-
  call bytes2human(nbytes_allocated, sizestring)
  if (nprocs == 1) then
     write(iprint, '(a)') '------------------------------------------------------------'
@@ -114,7 +111,6 @@ subroutine deallocate_memory(part_only)
  use photoevap, only:deallocate_photoevap
 #endif
  use allocutils, only:nbytes_allocated
- use omp_cache, only:deallocate_cache
 
  logical, optional, intent(in) :: part_only
  logical :: part_only_
@@ -133,8 +129,6 @@ subroutine deallocate_memory(part_only)
 #endif
  endif
 
- call deallocate_cache
-
  nbytes_allocated = 0
  call update_max_sizes(0)
 

diff --git a/src/main/omp_cache.f90 b/src/main/omp_cache.f90