Skip to content

Commit

Permalink
Merge branch 'master' of github.com:danieljprice/phantom
Browse files Browse the repository at this point in the history
  • Loading branch information
@danieljprice
danieljprice committed Aug 14, 2020
2 parents 1829b02 + b9c5d57 commit 1afdd3c
Show file tree
Hide file tree
Showing 21 changed files with 719 additions and 234 deletions.
13 changes: 12 additions & 1 deletion build/Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
6 changes: 4 additions & 2 deletions src/main/inject_asteroidwind.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
9 changes: 7 additions & 2 deletions src/main/part.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -592,6 +592,11 @@ subroutine init_part
dustgasprop(:,:) = 0.
VrelVf(:) = 0.
#endif
#ifdef IND_TIMESTEPS
ibin(:) = 0
ibin_old(:) = 0
ibin_wake(:) = 0
#endif

end subroutine init_part

Expand Down Expand Up @@ -806,8 +811,8 @@ pure integer(kind=1) function isetphase(itype,iactive)
integer, intent(in) :: itype
logical, intent(in) :: iactive

if ((set_boundaries_to_active .and. itype==iboundary) .or. &
(iactive .and. itype/=iboundary) ) then
if ((set_boundaries_to_active .and. iamboundary(itype)) .or. &
(iactive .and. .not.iamboundary(itype)) ) then
isetphase = int(itype,kind=1)
else
isetphase = -int(abs(itype),kind=1)
Expand Down
266 changes: 142 additions & 124 deletions src/main/readwrite_dumps.F90
View file
Original file line number Diff line number Diff line change
@@ -1,166 +1,184 @@
!--------------------------------------------------------------------------!
! 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 readwrite_dumps

use readwrite_dumps_fortran
!
! readwrite_dumps
!
! :References: None
!
! :Owner: Daniel Price
!
! :Runtime parameters: None
!
! :Dependencies: boundary, io, readwrite_dumps_fortran,
! readwrite_dumps_hdf5
!

use readwrite_dumps_fortran
#ifdef HDF5
use readwrite_dumps_hdf5
use readwrite_dumps_hdf5
#endif

implicit none
implicit none

public :: write_smalldump,write_fulldump,read_smalldump,read_dump,write_gadgetdump
public :: write_smalldump,write_fulldump,read_smalldump,read_dump,write_gadgetdump

logical, pointer, public :: opened_full_dump ! for use in analysis files if user wishes to skip small dumps
logical, pointer, public :: dt_read_in ! to determine if dt has been read in so that ibin & ibinold can be set on restarts
logical, pointer, public :: opened_full_dump ! for use in analysis files if user wishes to skip small dumps
logical, pointer, public :: dt_read_in ! to determine if dt has been read in so that ibin & ibinold can be set on restarts

integer, parameter, public :: is_small_dump = 1978
integer, parameter, public :: is_not_mhd = 1979
integer, parameter, public :: is_small_dump = 1978
integer, parameter, public :: is_not_mhd = 1979

procedure(read_dump_fortran), pointer :: read_dump => null()
procedure(read_smalldump_fortran), pointer :: read_smalldump => null()
procedure(write_smalldump_fortran), pointer :: write_smalldump => null()
procedure(write_fulldump_fortran), pointer :: write_fulldump => null()
procedure(read_dump_fortran), pointer :: read_dump => null()
procedure(read_smalldump_fortran), pointer :: read_smalldump => null()
procedure(write_smalldump_fortran), pointer :: write_smalldump => null()
procedure(write_fulldump_fortran), pointer :: write_fulldump => null()


contains

subroutine init_readwrite_dumps()
subroutine init_readwrite_dumps()

logical :: lhdf5
logical :: lhdf5

#ifdef HDF5
lhdf5 = .true. ! we can potentially use both file formats if needed

if (lhdf5) then
read_dump => read_dump_hdf5
read_smalldump => read_smalldump_hdf5
write_smalldump => write_smalldump_hdf5
write_fulldump => write_fulldump_hdf5

opened_full_dump => opened_full_dump_hdf5
dt_read_in => dt_read_in_hdf5
else
read_dump => read_dump_fortran
read_smalldump => read_smalldump_fortran
write_smalldump => write_smalldump_fortran
write_fulldump => write_fulldump_fortran

opened_full_dump => opened_full_dump_fortran
dt_read_in => dt_read_in_fortran
endif
#else
lhdf5 = .false.
lhdf5 = .true. ! we can potentially use both file formats if needed

if (lhdf5) then
read_dump => read_dump_hdf5
read_smalldump => read_smalldump_hdf5
write_smalldump => write_smalldump_hdf5
write_fulldump => write_fulldump_hdf5

opened_full_dump => opened_full_dump_hdf5
dt_read_in => dt_read_in_hdf5
else
read_dump => read_dump_fortran
read_smalldump => read_smalldump_fortran
write_smalldump => write_smalldump_fortran
write_fulldump => write_fulldump_fortran

opened_full_dump => opened_full_dump_fortran
dt_read_in => dt_read_in_fortran
endif
#else
lhdf5 = .false.

read_dump => read_dump_fortran
read_smalldump => read_smalldump_fortran
write_smalldump => write_smalldump_fortran
write_fulldump => write_fulldump_fortran

opened_full_dump => opened_full_dump_fortran
dt_read_in => dt_read_in_fortran
#endif

return
return

end subroutine init_readwrite_dumps
end subroutine init_readwrite_dumps

!--------------------------------------------------------------------
!+
! subroutine to write output to full dump file
! in GADGET format
!+
!-------------------------------------------------------------------
subroutine write_gadgetdump(dumpfile,t,xyzh,particlemass,vxyzu,rho,utherm,npart)
use io, only:iprint,idump,real4
!--------------------------------------------------------------------
!+
! subroutine to write output to full dump file
! in GADGET format
!+
!-------------------------------------------------------------------
subroutine write_gadgetdump(dumpfile,t,xyzh,particlemass,vxyzu,rho,utherm,npart)
use io, only:iprint,idump,real4
#ifdef PERIODIC
use boundary, only:dxbound
use boundary, only:dxbound
#endif
real, intent(in) :: t,particlemass,utherm
character(len=*), intent(in) :: dumpfile
integer, intent(in) :: npart
real, intent(in) :: xyzh(:,:),vxyzu(:,:)
real, intent(in) :: rho(:)

integer(kind=4) :: particleid(size(rho))
integer :: npartoftype(6),nall(6),ncrap(6)
real(kind=8) :: massoftype(6)
real(kind=8) :: time,boxsize
real(kind=8), parameter :: dumz = 0.d0
real(kind=4) :: unused(15)
integer, parameter :: iflagsfr = 0, iflagfeedback = 0, iflagcool = 0
integer, parameter :: nfiles = 1
integer :: ierr,i,j
!
!--open dumpfile
!
write(iprint,"(/,/,'--------> TIME = ',f12.4,"// &
real, intent(in) :: t,particlemass,utherm
character(len=*), intent(in) :: dumpfile
integer, intent(in) :: npart
real, intent(in) :: xyzh(:,:),vxyzu(:,:)
real, intent(in) :: rho(:)

integer(kind=4) :: particleid(size(rho))
integer :: npartoftype(6),nall(6),ncrap(6)
real(kind=8) :: massoftype(6)
real(kind=8) :: time,boxsize
real(kind=8), parameter :: dumz = 0.d0
real(kind=4) :: unused(15)
integer, parameter :: iflagsfr = 0, iflagfeedback = 0, iflagcool = 0
integer, parameter :: nfiles = 1
integer :: ierr,i,j
!
!--open dumpfile
!
write(iprint,"(/,/,'--------> TIME = ',f12.4,"// &
"': full dump written to file ',a,' <--------',/)") t,trim(dumpfile)

write(iprint,*) 'writing to unit ',idump
open(unit=idump,file=dumpfile,status='replace',form='unformatted',iostat=ierr)
if (ierr /= 0) then
write(iprint,*) 'error: can''t create new dumpfile ',trim(dumpfile)
stop
endif

npartoftype(:) = 0
npartoftype(1) = npart
nall(:) = npartoftype(:)
ncrap(:) = 0
time = t
write(iprint,*) 'writing to unit ',idump
open(unit=idump,file=dumpfile,status='replace',form='unformatted',iostat=ierr)
if (ierr /= 0) then
write(iprint,*) 'error: can''t create new dumpfile ',trim(dumpfile)
stop
endif

npartoftype(:) = 0
npartoftype(1) = npart
nall(:) = npartoftype(:)
ncrap(:) = 0
time = t
#ifdef PERIODIC
boxsize = dxbound
boxsize = dxbound
#else
boxsize = 0.
boxsize = 0.
#endif

massoftype(:) = 0.
massoftype(1) = particlemass
unused(:) = 0
massoftype(:) = 0.
massoftype(1) = particlemass
unused(:) = 0

do i=1,npart
particleid(i) = i
enddo
write(idump,iostat=ierr) npartoftype(1:6),massoftype(1:6),time,dumz, &
do i=1,npart
particleid(i) = i
enddo
write(idump,iostat=ierr) npartoftype(1:6),massoftype(1:6),time,dumz, &
iflagsfr,iflagfeedback,nall(1:6),iflagcool,nfiles,boxsize, &
dumz,dumz,dumz,iflagsfr,iflagsfr,ncrap(1:6),iflagsfr,unused(:)

write(idump,iostat=ierr) ((real4(xyzh(j,i)),j=1,3),i=1,npart)
if (ierr /= 0) then
print*,' error writing positions'
return
endif
write(idump,iostat=ierr) ((real4(vxyzu(j,i)),j=1,3),i=1,npart)
if (ierr /= 0) then
print*,' error writing velocities'
return
endif
write(idump,iostat=ierr) (particleid(i),i=1,npart)
if (ierr /= 0) then
print*,' error writing particle ID'
return
endif
if (size(vxyzu(:,1)) >= 4) then
write(idump,iostat=ierr) (real4(vxyzu(4,i)),i=1,npart)
else
write(idump,iostat=ierr) (real4(utherm),i=1,npart)
endif
if (ierr /= 0) then
print*,' error writing utherm'
return
endif
write(idump,iostat=ierr) (real4(rho(i)),i=1,npart)
if (ierr /= 0) then
print*,' error writing rho'
return
endif
write(idump,iostat=ierr) (real4(xyzh(4,i)),i=1,npart)
if (ierr /= 0) then
print*,' error writing h'
return
endif
print*,' finished writing file -- OK'

write(idump,iostat=ierr) ((real4(xyzh(j,i)),j=1,3),i=1,npart)
if (ierr /= 0) then
print*,' error writing positions'
return
endif
write(idump,iostat=ierr) ((real4(vxyzu(j,i)),j=1,3),i=1,npart)
if (ierr /= 0) then
print*,' error writing velocities'
return
endif
write(idump,iostat=ierr) (particleid(i),i=1,npart)
if (ierr /= 0) then
print*,' error writing particle ID'
return
endif
if (size(vxyzu(:,1)) >= 4) then
write(idump,iostat=ierr) (real4(vxyzu(4,i)),i=1,npart)
else
write(idump,iostat=ierr) (real4(utherm),i=1,npart)
endif
if (ierr /= 0) then
print*,' error writing utherm'
return
end subroutine write_gadgetdump
endif
write(idump,iostat=ierr) (real4(rho(i)),i=1,npart)
if (ierr /= 0) then
print*,' error writing rho'
return
endif
write(idump,iostat=ierr) (real4(xyzh(4,i)),i=1,npart)
if (ierr /= 0) then
print*,' error writing h'
return
endif
print*,' finished writing file -- OK'

return
end subroutine write_gadgetdump

end module readwrite_dumps
12 changes: 12 additions & 0 deletions src/main/readwrite_dumps_common.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,18 @@
! http://phantomsph.bitbucket.io/ !
!--------------------------------------------------------------------------!
module readwrite_dumps_common
!
! readwrite_dumps_common
!
! :References: None
!
! :Owner: Daniel Mentiplay
!
! :Runtime parameters: None
!
! :Dependencies: dim, dump_utils, eos, gitinfo, io, options, part,
! sphNGutils
!
use dump_utils, only:lenid
implicit none

Expand Down

0 comments on commit 1afdd3c

Please sign in to comment.