Merge pull request #41 from becnealon/master

Splitting/merging

build/Makefile

@@ -2109,12 +2109,23 @@ cleanstructutils: cleanstruct2struct cleanget_struct_slope
 #
 .PHONY: splitpart
 splitpart:
-	${MAKE} moddump MODFILE="splitpart.f90 moddump_splitpart.f90"\
+	${MAKE} moddump MODFILE="utils_splitmerge.f90 splitpart.f90 moddump_splitpart.f90"\
         MODDUMPBIN=$@
 
 cleansplitpart:
 	rm -f $(BINDIR)/splitpart
 
+#------------------------------------------------------
+# particle merging utility (this is a moddump compiled as a standalone utility)
+#
+.PHONY: mergepart
+mergepart:
+	${MAKE} moddump MODFILE="utils_getneighbours.F90 utils_splitmerge.f90 splitpart.f90 moddump_mergepart.f90"\
+        MODDUMPBIN=$@
+
+cleanmergepart:
+	rm -f $(BINDIR)/mergepart
+
 #----------------------------------------------------
 # utility to calculate divv from a dump file
 # compile using all phantom files

src/main/inject_asteroidwind.f90

@@ -77,7 +77,7 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,&
  integer :: i,ipart,npinject,seed,pt,test_integral
  real    :: dmdt,dndt,rasteroid,h,u,speed,inject_this_step
  real    :: m1,m2,mu,period,r,q
- real    :: phi,theta,mod_time,dt,func_now
+ real    :: phi,theta,mod_time,dt,func_now, phi_facing
 
  real, save :: have_injected,func_old,t_old
  real, save :: semia, ra, rp, ecc
@@ -173,9 +173,11 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,&
 
 !
 !-- Randomly inject particles around the asteroids outer 'radius'
+!-- Only inject them on the side that is facing the central sink
 !
  do i=1,npinject
-    phi       = ran2(seed)*twopi
+    phi_facing = atan2(xyzmh_ptmass(2,pt),xyzmh_ptmass(1,pt)) + 0.5*pi
+    phi       = phi_facing + ran2(seed)*pi
     theta     = ran2(seed)*pi
     xyz       = r2 + (/rasteroid*cos(phi)*sin(theta),rasteroid*sin(phi)*sin(theta),rasteroid*cos(theta)/)
     vxyz      = (1.-vlag/100)*speed*vhat

src/utils/analysis_clumpfind.F90

@@ -18,11 +18,11 @@ module analysis
 !
 ! :Runtime parameters: None
 !
-! :Dependencies: boundary, dim, getneigbours, part, prompting, ptmass,
+! :Dependencies: boundary, dim, getneighbours, part, prompting, ptmass,
 !   readwrite_dumps, sortutils, units
 !
  use dim,             only:maxp
- use getneigbours,    only:generate_neighbour_lists, read_neighbours, write_neighbours, &
+ use getneighbours,    only:generate_neighbour_lists, read_neighbours, write_neighbours, &
                            neighcount,neighb,neighmax
  implicit none
  character(len=20), parameter, public :: analysistype = 'clumpfind'

src/utils/analysis_disc_stresses.f90

@@ -14,10 +14,10 @@ module analysis
 !
 ! :Runtime parameters: None
 !
-! :Dependencies: dim, eos, getneigbours, io, kernel, part, physcon,
+! :Dependencies: dim, eos, getneighbours, io, kernel, part, physcon,
 !   prompting, units
 !
- use getneigbours,    only:generate_neighbour_lists, read_neighbours, write_neighbours, &
+ use getneighbours,    only:generate_neighbour_lists, read_neighbours, write_neighbours, &
                            neighcount,neighb,neighmax
  implicit none
  character(len=20), parameter, public :: analysistype = 'disc_stresses'

src/utils/analysis_getneighbours.F90

@@ -14,9 +14,9 @@ module analysis
 !
 ! :Runtime parameters: None
 !
-! :Dependencies: getneigbours
+! :Dependencies: getneighbours
 !
- use getneigbours,    only:generate_neighbour_lists, read_neighbours, write_neighbours, &
+ use getneighbours,    only:generate_neighbour_lists, read_neighbours, write_neighbours, &
                            neighcount,neighb,neighmax
  implicit none
  character(len=20), parameter, public :: analysistype = 'getneighbours'

src/utils/analysis_velocitydispersion_vs_scale.f90

@@ -14,10 +14,10 @@ module analysis
 !
 ! :Runtime parameters: None
 !
-! :Dependencies: dim, eos, getneigbours, io, part, physcon, prompting,
+! :Dependencies: dim, eos, getneighbours, io, part, physcon, prompting,
 !   sortutils
 !
- use getneigbours,    only:generate_neighbour_lists, read_neighbours, write_neighbours, &
+ use getneighbours,    only:generate_neighbour_lists, read_neighbours, write_neighbours, &
                            neighcount,neighb,neighmax,meanneigh
  implicit none
  character(len=27), parameter, public :: analysistype = 'velocitydispersion_vs_scale'

src/utils/analysis_velocityshear.F90

@@ -19,9 +19,9 @@ module analysis
 !
 ! :Runtime parameters: None
 !
-! :Dependencies: dim, getneigbours, kernel, part
+! :Dependencies: dim, getneighbours, kernel, part
 !
- use getneigbours,    only:generate_neighbour_lists, read_neighbours, write_neighbours, &
+ use getneighbours,    only:generate_neighbour_lists, read_neighbours, write_neighbours, &
                            neighcount,neighb,neighmax
  implicit none
  character(len=20), parameter, public :: analysistype = 'velocityshear'

src/utils/moddump_mergepart.f90

@@ -0,0 +1,64 @@
+!--------------------------------------------------------------------------!
+! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. !
+! Copyright (c) 2007-2020 The Authors (see AUTHORS)                        !
+! See LICENCE file for usage and distribution conditions                   !
+! http://phantomsph.bitbucket.io/                                          !
+!--------------------------------------------------------------------------!
+!+
+!  MODULE: moddump
+!
+!  DESCRIPTION:
+!  merges particles; input simulation with npart, get npart/nchild back
+!
+!  REFERENCES: Vacondio et al. 2013
+!
+!  OWNER: Rebecca Nealon
+!
+!  $Id$
+!
+!  RUNTIME PARAMETERS: None
+!
+!  DEPENDENCIES: io, splitpart
+!+
+!--------------------------------------------------------------------------
+module moddump
+ implicit none
+ integer, parameter :: nchild = 2
+
+contains
+
+subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
+ use splitpart, only:merge_all_particles
+ use part,      only:igas,kill_particle,delete_dead_or_accreted_particles
+ use part,      only:isdead_or_accreted,copy_particle
+ use io,        only:fatal,error
+ integer, intent(inout) :: npart
+ integer, intent(inout) :: npartoftype(:)
+ real,    intent(inout) :: massoftype(:)
+ real,    intent(inout) :: xyzh(:,:),vxyzu(:,:)
+ integer :: i,nactive
+
+ !-- how many active particles
+ nactive = 0
+ do i = 1,npart
+   if (.not.isdead_or_accreted(xyzh(4,i))) then
+     nactive = nactive + 1
+   else
+     call kill_particle(i,npartoftype)
+   endif
+ enddo
+
+ if (nactive < npart) then
+   call delete_dead_or_accreted_particles(npart,npartoftype)
+   print*,' discarding inactive particles'
+ endif
+
+ ! Merge 'em!
+ call merge_all_particles(npart,npartoftype,massoftype,xyzh,vxyzu, &
+                                nchild,nactive)
+
+ print*,' new npart = ',npart
+
+end subroutine modify_dump
+
+end module moddump

src/utils/moddump_splitpart.f90

@@ -14,25 +14,41 @@ module moddump
 !
 ! :Runtime parameters: None
 !
-! :Dependencies: io, splitpart
+! :Dependencies: io, splitpart, injectutils
 !
  implicit none
- integer, parameter :: nchildren = 8
+ integer            :: nchild = 12
+ integer, parameter :: lattice_type = 0 ! 0 for lattice, 1 for random
+ integer, parameter :: ires = 1         ! use 12 particles per sphere
 
 contains
 
 subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
- use splitpart, only:split_particles
- use io,        only:fatal
+ use splitpart,    only:split_all_particles
+ use io,           only:fatal,error
+ use injectutils,  only:get_parts_per_sphere
+ use part,         only:delete_dead_or_accreted_particles
  integer, intent(inout) :: npart
  integer, intent(inout) :: npartoftype(:)
  real,    intent(inout) :: massoftype(:)
  real,    intent(inout) :: xyzh(:,:),vxyzu(:,:)
  integer :: ierr
 
- call split_particles(nchildren,npart,npartoftype,xyzh,massoftype,ierr)
- if (ierr /= 0) call fatal('moddump','could not split particles')
- print*,' got npart = ',npart
+ ierr = 0
+
+ !-- if using the regular grid, set nchild to get desired resolution
+ if (lattice_type == 0) then
+   nchild = get_parts_per_sphere(ires) + 1
+ endif
+
+ !-- don't split accreted particles
+ call delete_dead_or_accreted_particles(npart,npartoftype)
+
+ ! Split 'em!
+ call split_all_particles(npart,npartoftype,massoftype,xyzh,vxyzu, &
+                                nchild,lattice_type,ires)
+
+ print*,' new npart = ',npart
 
 end subroutine modify_dump