Merge pull request danieljprice#620 from Yrisch/sinktree_matrix

Sinks in the kd tree

.github/workflows/mpi.yml

@@ -26,6 +26,7 @@ jobs:
           - ['testdust', 'dust']
           - ['testgr', 'gr']
           - ['testgrav', 'gravity ptmass']
+          - ['testsinktree','gravity ptmass']
           - ['testgrowth', 'dustgrowth']
           - ['testnimhd', 'nimhd']
           - ['test2', '']

.github/workflows/test.yml

@@ -33,6 +33,7 @@ jobs:
           - ['testdust', 'dust']
           - ['testgr', 'gr ptmass']
           - ['testgrav', 'gravity ptmass setstar']
+          - ['testsinktree','gravity ptmass']
           - ['testgrowth', 'dustgrowth']
           - ['testnimhd', 'nimhd']
           - ['test2', '']

build/Makefile

@@ -283,6 +283,10 @@ ifdef SRCTURB
     FPPFLAGS += -DDRIVING
 endif
 
+ifeq ($(SINKTREE), yes)
+    FPPFLAGS += -DSINKTREE
+endif
+
 ifeq ($(KROME), krome)
     FPPFLAGS += -DKROME
 ifeq ($(SYSTEM), ifort)

build/Makefile_setups

@@ -1008,6 +1008,15 @@ ifeq ($(SETUP), testgrav)
     KNOWN_SETUP=yes
 endif
 
+ifeq ($(SETUP), testsinktree)
+#   self-gravity unit tests
+    GRAVITY=yes
+    SINKTREE=yes
+    CONST_ARTRES=yes
+    CURLV=yes
+    KNOWN_SETUP=yes
+endif
+
 ifeq ($(SETUP), testdust)
 #   dust unit tests
     PERIODIC=yes

src/main/H2regions.f90

@@ -187,7 +187,7 @@ subroutine HII_feedback(nptmass,npart,xyzh,xyzmh_ptmass,vxyzu,isionised,dt)
  use sortutils,  only:Knnfunc,set_r2func_origin,r2func_origin
  use physcon,    only:pc,pi
  use timing,     only:get_timings,increment_timer,itimer_HII
- use dim,        only:maxvxyzu
+ use dim,        only:maxvxyzu,maxpsph
  use units,      only:utime
  integer,          intent(in)    :: nptmass,npart
  real,             intent(in)    :: xyzh(:,:)
@@ -240,6 +240,7 @@ subroutine HII_feedback(nptmass,npart,xyzh,xyzmh_ptmass,vxyzu,isionised,dt)
        enddo
        do k=1,nneigh
           j = listneigh(k)
+          if (j > maxpsph) cycle
           if (.not. isdead_or_accreted(xyzh(4,j))) then
              ! ionising photons needed to fully ionise the current particle
              DNdot = log10((((pmass*ar*rhoh(xyzh(4,j),pmass))/(mH**2))/utime))

src/main/config.F90

@@ -128,7 +128,7 @@ module dim
                                    radensumden
 
  ! fsum
- integer, parameter :: fsumvars = 23 ! Number of scalars in fsum
+ integer, parameter :: fsumvars = 25 ! Number of scalars in fsum
  integer, parameter :: fsumarrs = 5  ! Number of arrays  in fsum
  integer, parameter :: maxfsum  = fsumvars + &                  ! Total number of values
                                   fsumarrs*(maxdusttypes-1) + &
@@ -137,7 +137,7 @@ module dim
 ! xpartveci
  integer, parameter :: maxxpartvecidens = 14 + radenxpartvetden
 
- integer, parameter :: maxxpartvecvars = 62 ! Number of scalars in xpartvec
+ integer, parameter :: maxxpartvecvars = 63 ! Number of scalars in xpartvec
  integer, parameter :: maxxpartvecarrs = 2  ! Number of arrays in xpartvec
  integer, parameter :: maxxpartvecGR   = 33 ! Number of GR values in xpartvec (1 for dens, 16 for gcov, 16 for gcon)
  integer, parameter :: maxxpartveciforce = maxxpartvecvars + &              ! Total number of values
@@ -185,7 +185,7 @@ module dim
 
  ! Maximum number of particle types
  !
- integer, parameter :: maxtypes = 7 + 2*maxdustlarge - 1
+ integer, parameter :: maxtypes = 8 + 2*maxdustlarge - 1
 
  !
  ! Number of dimensions, where it is needed
@@ -339,6 +339,16 @@ module dim
 #endif
  integer :: maxp_apr = 0
 
+!--------------------
+! Sink in tree methods
+!--------------------
+#ifdef SINKTREE
+ logical, parameter :: use_sinktree = .true.
+#else
+ logical, parameter :: use_sinktree = .false.
+#endif
+ integer :: maxpsph = 0
+
 !--------------------
 ! individual timesteps
 !--------------------
@@ -380,6 +390,11 @@ subroutine update_max_sizes(n,ntot)
     maxp_apr = maxp
  endif
 
+ maxpsph = maxp
+ if (use_sinktree) then
+    maxp = n+maxptmass
+ endif
+
  if (use_krome) maxp_krome = maxp
 
  if (h2chemistry) maxp_h2 = maxp

src/main/dens.F90

@@ -605,7 +605,7 @@ 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
  use part,     only:massoftype,iradxi,aprmassoftype
- use dim,      only:ndivcurlv,gravity,maxp,nalpha,use_dust,do_radiation,use_apr
+ use dim,      only:ndivcurlv,gravity,maxp,nalpha,use_dust,do_radiation,use_apr,maxpsph
  use options,  only:implicit_radiation
  integer,      intent(in)    :: i
  real,         intent(in)    :: xpartveci(:)
@@ -668,6 +668,7 @@ pure subroutine get_density_sums(i,xpartveci,hi,hi1,hi21,iamtypei,iamgasi,iamdus
     j = listneigh(n)
     !--do self contribution separately to avoid problems with 1/sqrt(0.)
     if ((ignoreself) .and. (j==i)) cycle loop_over_neigh
+    if (j > maxpsph) cycle loop_over_neigh
 
     if (ifilledneighcache .and. n <= isizeneighcache) then
        rij2 = dxcache(1,n)
@@ -1303,7 +1304,7 @@ end subroutine compute_hmax
 !--------------------------------------------------------------------------
 subroutine start_cell(cell,iphase,xyzh,vxyzu,fxyzu,fext,Bevol,rad,apr_level)
  use io,          only:fatal
- use dim,         only:maxp,maxvxyzu,do_radiation,use_apr
+ use dim,         only:maxp,maxvxyzu,do_radiation,use_apr,maxpsph
  use part,        only:maxphase,get_partinfo,mhd,igas,iamgas,&
                        iamboundary,ibasetype,iradxi
 
@@ -1325,7 +1326,7 @@ subroutine start_cell(cell,iphase,xyzh,vxyzu,fxyzu,fext,Bevol,rad,apr_level)
  over_parts: do ip = inoderange(1,cell%icell),inoderange(2,cell%icell)
     i = inodeparts(ip)
 
-    if (i < 0) then
+    if (i < 0 .or. i > maxpsph) then
        cycle over_parts
     endif
 

src/main/deriv.F90

@@ -224,7 +224,7 @@ end subroutine derivs
 !  does not work for sink GR yet
 !+
 !--------------------------------------
-subroutine get_derivs_global(tused,dt_new,dt)
+subroutine get_derivs_global(tused,dt_new,dt,icall)
  use part,         only:npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,&
                         Bevol,dBevol,rad,drad,radprop,dustprop,ddustprop,filfac,&
                         dustfrac,ddustevol,eos_vars,pxyzu,dens,metrics,dustevol,gr,&
@@ -235,14 +235,18 @@ subroutine get_derivs_global(tused,dt_new,dt)
  use metric_tools, only:init_metric
  real(kind=4), intent(out), optional :: tused
  real,         intent(out), optional :: dt_new
- real,         intent(in), optional  :: dt  ! optional argument needed to test implicit radiation routine
+ real,         intent(in),  optional :: dt  ! optional argument needed to test implicit radiation routine
+ integer     , intent(in),  optional :: icall
  real(kind=4) :: t1,t2
- real :: dtnew
- real :: time,dti
+ real    :: dtnew
+ real    :: time,dti
+ integer :: icalli
 
  time = 0.
  dti = 0.
+ icalli = 1
  if (present(dt)) dti = dt
+ if (present(icall)) icalli = icall
  call getused(t1)
  ! update conserved quantities in the GR code
  if (gr) then
@@ -251,7 +255,7 @@ subroutine get_derivs_global(tused,dt_new,dt)
  endif
 
  ! evaluate derivatives
- call derivs(1,npart,npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,Bevol,dBevol,&
+ call derivs(icalli,npart,npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,Bevol,dBevol,&
              rad,drad,radprop,dustprop,ddustprop,dustevol,ddustevol,filfac,dustfrac,&
              eos_vars,time,dti,dtnew,pxyzu,dens,metrics,apr_level)
  call getused(t2)

src/main/energies.F90

@@ -65,7 +65,7 @@ module energies
 subroutine compute_energies(t)
  use dim,            only:maxp,maxvxyzu,maxalpha,maxtypes,mhd_nonideal,maxp_hard,&
                           lightcurve,use_dust,maxdusttypes,do_radiation,gr,use_krome,&
-                          use_apr
+                          use_apr,use_sinktree,maxpsph
  use part,           only:rhoh,xyzh,vxyzu,massoftype,npart,maxphase,iphase,&
                           alphaind,Bevol,divcurlB,iamtype,igamma,&
                           igas,idust,iboundary,istar,idarkmatter,ibulge,&
@@ -171,7 +171,7 @@ subroutine compute_energies(t)
  call initialise_ev_data(ev_data)
 
 !$omp parallel default(none) &
-!$omp shared(maxp,maxphase,maxalpha) &
+!$omp shared(maxp,maxphase,maxalpha,maxpsph) &
 !$omp shared(xyzh,vxyzu,pxyzu,rad,iexternalforce,npart,t,id) &
 !$omp shared(alphaind,massoftype,irealvisc,iu,aprmassoftype) &
 !$omp shared(ieos,gamma,nptmass,xyzmh_ptmass,vxyz_ptmass,xyzcom) &
@@ -351,7 +351,7 @@ subroutine compute_energies(t)
           call externalforce_vdependent(iexternalforce,xyzh(1:3,i),vxyzu(1:3,i),fdum,epottmpi)
           epoti = pmassi*epottmpi
        endif
-       if (nptmass > 0) then
+       if (nptmass > 0 .and. .not.use_sinktree) then ! No need to compute if sink in tree
           dumx = 0.
           dumy = 0.
           dumz = 0.
@@ -595,6 +595,7 @@ subroutine compute_energies(t)
           angy   = angy + pmassi*(zi*vxi - xi*vzi)
           angz   = angz + pmassi*(xi*vyi - yi*vxi)
 
+          if (use_sinktree) epot = epot + poten(i+maxpsph)
           angx   = angx + xyzmh_ptmass(ispinx,i)
           angy   = angy + xyzmh_ptmass(ispiny,i)
           angz   = angz + xyzmh_ptmass(ispinz,i)

src/main/evolve.F90

@@ -93,7 +93,8 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag)
                             xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dptmass,gravity,iboundary, &
                             fxyz_ptmass_sinksink,ntot,poten,ndustsmall,&
                             accrete_particles_outside_sphere,apr_level,aprmassoftype,&
-                            sf_ptmass,isionised,dsdt_ptmass,isdead_or_accreted
+                            sf_ptmass,isionised,dsdt_ptmass,isdead_or_accreted,&
+                            fxyz_ptmass_tree
  use part,             only:n_group,n_ingroup,n_sing,group_info,bin_info,nmatrix
  use quitdump,         only:quit
  use ptmass,           only:icreate_sinks,ptmass_create,ipart_rhomax,pt_write_sinkev,calculate_mdot, &
@@ -336,8 +337,8 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag)
                               new_ptmass=.true.,dtext=dtextforce)
        endif
        call get_force(nptmass,npart,0,1,time,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass,vxyz_ptmass,&
-                      fxyz_ptmass,dsdt_ptmass,0.,0.,dummy,.false.,sf_ptmass,bin_info,group_info,&
-                      nmatrix)
+                      fxyz_ptmass,fxyz_ptmass_tree,dsdt_ptmass,0.,0.,dummy,.false.,sf_ptmass,&
+                      bin_info,group_info,nmatrix)
        if (ipart_createseeds /= 0) ipart_createseeds = 0 ! reset pointer to zero
        if (ipart_createstars /= 0) ipart_createstars = 0 ! reset pointer to zero
        dummy = 0