Fmm prec plot

src/main/kdtree.F90

@@ -56,7 +56,6 @@ module kdtree
  public :: maketreeglobal
  public :: empty_tree
  public :: compute_M2L,expand_fgrav_in_taylor_series
-
  integer, public :: maxlevel_indexed, maxlevel
 
  type kdbuildstack

src/main/neigh_kdtree.f90

@@ -36,7 +36,7 @@ module neighkdtree
  integer(kind=8), public            :: ncells
  real, public                       :: dxcell
  real, public                       :: dcellx = 0.,dcelly = 0.,dcellz = 0.
- logical, public                    :: force_dual_walk
+ logical, public                    :: use_dualtree = .true.
  integer                            :: globallevel,refinelevels
 
  public :: allocate_neigh, deallocate_neigh
@@ -279,12 +279,12 @@ subroutine get_neighbour_list(inode,mylistneigh,nneigh,xyzh,xyzcache,ixyzcachesi
  endif
 
  ! Find neighbours of this cell on this node
- if ((get_f .or. force_dual_walk) .and. (.not.mpi)) then
+ if (get_f .and. .not.(mpi) .and. use_dualtree) then
     call getneigh_dual(node,xpos,xsizei,rcuti,mylistneigh,nneigh,xyzcache,ixyzcachesize,&
-                       leaf_is_active,get_j,get_f,fgrav,inode)
+                          leaf_is_active,get_j,get_f,fgrav,inode)
  else
     call getneigh(node,xpos,xsizei,rcuti,mylistneigh,nneigh,xyzcache,ixyzcachesize,&
-                  leaf_is_active,get_j,get_f,fgrav)
+                     leaf_is_active,get_j,get_f,fgrav)
  endif
 
  if (get_f) f = fgrav + fgrav_global
@@ -326,9 +326,7 @@ subroutine write_options_tree(iunit)
  use part,         only:gravity
  integer, intent(in) :: iunit
 
- if (gravity) then
-    call write_inopt(tree_accuracy,'tree_accuracy','tree opening criterion (0.0-1.0)',iunit)
- endif
+ if (gravity) call write_inopt(tree_accuracy,'tree_accuracy','tree opening criterion (0.0-1.0)',iunit)
 
 end subroutine write_options_tree
 

src/main/writeheader.f90

@@ -158,7 +158,7 @@ subroutine write_header(icall,infile,evfile,logfile,dumpfile,ntot)
            6x,' Number of neighbours = ',i4,/)
 
     if (mhd)              write(iprint,"(1x,a)") 'Magnetic fields are ON, evolving B/rho with cleaning'
-    if (gravity)          write(iprint,"(1x,a)") 'Self-gravity is ON'
+    if (gravity)          write(iprint,"(1x,a)") 'Self-gravity is ON, please cite Bernard et al. (2026) and Price & Monaghan (2007)'
     if (h2chemistry)      write(iprint,"(1x,a)") 'H2 Chemistry is ON'
     if (use_dust) then
        if (use_dustfrac) then

src/setup/set_sphere.f90

@@ -24,6 +24,8 @@ module spherical
 
  public  :: set_sphere,set_ellipse,rho_func
 
+ integer,public :: iseed_mc = -1978
+
  integer, parameter :: &
    ierr_notinrange    = 1, &
    ierr_not_converged = 2, &
@@ -162,7 +164,7 @@ subroutine set_sphere_mc(id,master,rmin,rmax,hfact,np_requested,np,xyzh, &
  integer(kind=8),  intent(inout) :: nptot
  logical,          intent(in)    :: verbose
  procedure(mask_prototype) :: mask
- integer :: i,npin,iseed,maxp
+ integer :: i,npin,maxp
  real    :: vol_sphere,rr,phi,theta,mr,dir(3)
  real    :: sintheta,costheta,sinphi,cosphi,psep
  integer(kind=8) :: iparttot
@@ -172,24 +174,23 @@ subroutine set_sphere_mc(id,master,rmin,rmax,hfact,np_requested,np,xyzh, &
  vol_sphere = 4./3.*pi*(rmax**3 - rmin**3)
  ! use mean particle spacing to set initial smoothing lengths
  psep = (vol_sphere/real(np_requested))**(1./3.)
- iseed = -1978
  maxp  = size(xyzh(1,:))
  ierr  = 1
 
  do i=npin+1,npin+np_requested,2
     !
     ! get random mass coordinate i.e. m(r)
     !
-    mr = ran2(iseed)
+    mr = ran2(iseed_mc)
     !
     ! invert to get mass coordinate from radial coordinate, i.e. r(m)
     !
     rr = (rmin**3+mr*(rmax**3-rmin**3))**(1./3.) ! uniform density
     !
     ! get a random position on sphere
     !
-    phi = 2.*pi*(ran2(iseed) - 0.5)
-    costheta = 2.*ran2(iseed) - 1.
+    phi = 2.*pi*(ran2(iseed_mc) - 0.5)
+    costheta = 2.*ran2(iseed_mc) - 1.
     theta    = acos(costheta)
     sintheta = sin(theta)
     sinphi   = sin(phi)