Merge branch 'master' of https://bitbucket.org/danielprice/phantom

.mailmap

@@ -20,12 +20,13 @@ Terrence Tricco <ttricco@cita.utoronto.ca>
 <ttricco@cita.utoronto.ca> <tricco@rin.(none)>
 <ttricco@cita.utoronto.ca> <ttricco@g2.hpc.swin.edu.au>
 <ttricco@cita.utoronto.ca> <terrence.tricco@monash.edu>
-James Wurster <j.wurster@exeter.ac.uk> jameswurster <jameswurster@bitbucket.org>
-James Wurster <j.wurster@exeter.ac.uk> jameswurster <james.wurster@monash.edu>
-James Wurster <j.wurster@exeter.ac.uk> <jwurster@g2.hpc.swin.edu.au>
-James Wurster <j.wurster@exeter.ac.uk> James Wurster <james.wurster@monash.edu>
-James Wurster <j.wurster@exeter.ac.uk> James <james.wurster@monash.edu>
-James Wurster <j.wurster@exeter.ac.uk> James Wurster <jw754@service2.ice.astro.ex.ac.uk>
+James Wurster <jhw5@st-andrews.ac.uk> James Wurster <j.wurster@exeter.ac.uk>
+James Wurster <jhw5@st-andrews.ac.uk> jameswurster <jameswurster@bitbucket.org>
+James Wurster <jhw5@st-andrews.ac.uk> jameswurster <james.wurster@monash.edu>
+James Wurster <jhw5@st-andrews.ac.uk> <jwurster@g2.hpc.swin.edu.au>
+James Wurster <jhw5@st-andrews.ac.uk> James Wurster <james.wurster@monash.edu>
+James Wurster <jhw5@st-andrews.ac.uk> James <james.wurster@monash.edu>
+James Wurster <jhw5@st-andrews.ac.uk> James Wurster <jw754@service2.ice.astro.ex.ac.uk>
 Ben Ayliffe <b.a.ayliffe@gmail.com>
 Owen Kaluza <owen.kaluza@monash.edu> okaluza <okaluza@msgln5.its.monash.edu.au>
 Owen Kaluza <owen.kaluza@monash.edu> okaluza <okaluza@msgln6.its.monash.edu.au>

AUTHORS

@@ -7,7 +7,7 @@
 #-------------------------------------------------------#
 Daniel Price <daniel.price@monash.edu>
 Conrad Chan <conrad.chan@monash.edu>
-James Wurster <j.wurster@exeter.ac.uk>
+James Wurster <jhw5@st-andrews.ac.uk>
 Daniel Mentiplay <daniel.mentiplay@monash.edu>
 Arnaud Vericel <arnaud.vericel@univ-lyon1.fr>
 Mark Hutchison <markahutch@gmail.com>
@@ -36,6 +36,7 @@ Nicolas Cuello <cuellonicolas@gmail.com>
 Sergei Biriukov <svbiriukov@gmail.com>
 Giulia Ballabio <giulia.ballabio2@studenti.unimi.it>
 Hauke Worpel <hauke@ursa.aip.de>
+Joe Fisher <jwfis1@student.monash.edu>
 Martina Toscani <mtoscani94@gmail.com>
 Orsola De Marco <orsola.demarco@mq.edu.au>
 Rebecca Nealon <rebecca.nealon@leicester.ac.uk>

docs/dustgrowth.rst

@@ -6,7 +6,7 @@ al. (in prep) <https://media.giphy.com/media/XIqCQx02E1U9W/giphy.>`__.
 
 .. important::
 
- Dust growth is not restrained to single star simulations!
+ Dust growth is not restricted to single star simulations!
  You can use dustgrowth around binaries, you can put planets in your
  disc, you can tilt and warp everything if you want and many more. Be
  creative :).
@@ -51,8 +51,8 @@ adds an option block after the dust section:
                   rsnow =        100.    ! snow line position in AU
                   Tsnow =        150.    ! snow line condensation temperature in K
                   vfrag =         15.    ! uniform fragmentation threshold in m/s
-                vfragin =       5.000    ! inward fragmentation threshold in m/s
-               vfragout =         15.    ! inward fragmentation threshold in m/s
+                vfragin =       5.000    ! internal fragmentation threshold in m/s
+               vfragout =         15.    ! external fragmentation threshold in m/s
            grainsizemin =       0.005    ! minimum allowed grain size in cm
 
 Here’s a brief description of each of them (remember that they are
@@ -97,7 +97,7 @@ occurs (if ifrag > 0 and isnow = 0)
                 vfragin =       5.000    ! inward fragmentation threshold in m/s
                vfragout =         15.    ! inward fragmentation threshold in m/s
 
-inner and outer fragmentation threshold velocities above which
+internal and external fragmentation threshold velocities above which
 fragmentation occurs (if ifrag > 0 and isnow > 0)
 
 ::

docs/examples.rst

@@ -12,3 +12,4 @@ This section contains some examples of physical problems that you can solve with
    star
    dustsettle
    dustgrowth
+   density

docs/hdf5.rst

@@ -310,9 +310,6 @@ Then you can access datasets like
 Plonk
 ~~~~~
 
-Plonk is a Python package for analysis and visualisation of Phantom
-data. It is open source and available at
+Plonk is a Python package for analysis and visualisation of SPH
+data. Plonk is open source and available at
 https://github.com/dmentipl/plonk.
-
-It uses compiled Fortran from Splash to do interpolation from the data
-to a pixel grid.

docs/index.rst

@@ -10,6 +10,10 @@ The Phantom SPH code, by Daniel Price.
 
 This code is designed to be an ultra-sleek, ultra-low-memory, code for high resolution SPH simulations.
 
+- Project homepage: http://phantomsph.bitbucket.io/
+- Code repository: https://bitbucket.org/danielprice/phantom/
+- Documentation: https://phantomsph.readthedocs.org/
+
 Contents
 --------
 

docs/plonk.rst

@@ -1,9 +1,21 @@
 Plonk
 =====
 
-Plonk is a Python package for analysis and visualisation of smoothed particle hydrodynamics data. It is open source.
-
-It uses Splash for visualisation.
+Plonk is a Python package for analysis and visualisation of smoothed particle
+hydrodynamics data. Plonk is open source.
 
 - Docs: https://plonk.readthedocs.io/
 - Repo: https://www.github.com/dmentipl/plonk
+
+Examples
+--------
+
+.. figure:: _static/images/plonk_render.png
+
+Density projection render of a snapshot in Cartesian coordinates, and polar
+coordinates. The data comes from a single Phantom dump.
+
+.. figure:: _static/images/plonk_accretion.png
+
+Mass accretion and accretion rate onto sink particles. The data comes from the
+Phantom sink `.ev` files.

src/main/force.F90

@@ -836,9 +836,9 @@ subroutine compute_forces(i,iamgasi,iamdusti,xpartveci,hi,hi1,hi21,hi41,gradhi,g
  real    :: dendissterm,dBdissterm,dudtresist,dpsiterm,pmassonrhoi
  real    :: gradpi,projsxi,projsyi,projszi
  real    :: gradp,projsx,projsy,projsz,Bxj,Byj,Bzj,Bj,Bj1,psij
- real    :: dpsitermj,grkernj,grgrkernj,autermj,avBtermj,vsigj,spsoundj
- real    :: gradpj,pro2j,projsxj,projsyj,projszj,sxxj,sxyj,sxzj,syyj,syzj,szzj,psitermj,dBrhoterm
- real    :: visctermisoj,visctermanisoj,enj,tempj,hj,mrhoj5,alphaj,pmassj,rho1j,vsigBj
+ real    :: grkernj,grgrkernj,autermj,avBtermj,vsigj,spsoundj
+ real    :: gradpj,pro2j,projsxj,projsyj,projszj,sxxj,sxyj,sxzj,syyj,syzj,szzj,dBrhoterm
+ real    :: visctermisoj,visctermanisoj,enj,tempj,hj,mrhoj5,alphaj,pmassj,rho1j
  real    :: rhoj,ponrhoj,prj,rhoav1
  real    :: hj1,hj21,q2j,qj,vwavej,divvj
  real    :: dvdxi(9),dvdxj(9)
@@ -990,9 +990,6 @@ subroutine compute_forces(i,iamgasi,iamdusti,xpartveci,hi,hi1,hi21,hi41,gradhi,g
  projsxj   = 0.
  projsyj   = 0.
  projszj   = 0.
- dpsitermj = 0.
- psitermj  = 0.
- vsigBj    = 0.
  dudtresist = 0.
  dpsiterm   = 0.
  fgravi = 0.

src/main/growth.F90

@@ -172,7 +172,7 @@ end subroutine print_growthinfo
 !+
 !-----------------------------------------------------------------------
 subroutine get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,dsdt)
- use part,            only:get_pmass,rhoh,idust,iamtype,iphase,isdead_or_accreted
+ use part,            only:rhoh,idust,iamtype,iphase,isdead_or_accreted,massoftype
  real, intent(in)     :: dustprop(:,:),dustgasprop(:,:)
  real, intent(in)     :: xyzh(:,:)
  real, intent(inout)  :: VrelVf(:),vxyzu(:,:)
@@ -192,7 +192,7 @@ subroutine get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,dsdt)
 
        if (iam == idust) then
 
-          rhod = rhoh(xyzh(4,i),get_pmass(idust,.false.))
+          rhod = rhoh(xyzh(4,i),massoftype(idust))
 
           call get_vrelonvfrag(xyzh(:,i),vxyzu(:,i),vrel,VrelVf(i),dustgasprop(:,i))
           !

src/main/readwrite_dumps_hdf5.F90

@@ -471,7 +471,7 @@ end subroutine write_dump
 !+
 !-------------------------------------------------------------------
 subroutine read_dump(dumpfile,tfile,hfactfile,idisk1,iprint,id,nprocs,ierr,headeronly,dustydisc)
- use boundary,       only:xmin,xmax,ymin,ymax,zmin,zmax
+ use boundary,       only:set_boundary
  use dim,            only:maxp,gravity,maxalpha,mhd,use_dust,use_dustgrowth, &
                           h2chemistry,store_temperature,nsinkproperties,maxp_hard
  use eos,            only:ieos,polyk,gamma,polyk2,qfacdisc,isink
@@ -483,7 +483,8 @@ subroutine read_dump(dumpfile,tfile,hfactfile,idisk1,iprint,id,nprocs,ierr,heade
                           nptmass,xyzmh_ptmass,vxyz_ptmass,ndustlarge,        &
                           ndustsmall,grainsize,graindens,Bextx,Bexty,Bextz,   &
                           alphaind,poten,Bxyz,Bevol,dustfrac,deltav,dustprop, &
-                          tstop,dustgasprop,VrelVf,temperature,abundance,ndusttypes
+                          tstop,dustgasprop,VrelVf,temperature,abundance,     &
+                          periodic,ndusttypes
 #ifdef IND_TIMESTEPS
  use part,           only:dt_in
 #endif
@@ -502,6 +503,7 @@ subroutine read_dump(dumpfile,tfile,hfactfile,idisk1,iprint,id,nprocs,ierr,heade
 
  real(kind=4), allocatable :: dtind(:)
 
+ real :: xmin,xmax,ymin,ymax,zmin,zmax
  character(len=200) :: fileident
  integer :: errors(5)
  logical :: smalldump,isothermal,ind_timesteps,const_av
@@ -603,6 +605,10 @@ subroutine read_dump(dumpfile,tfile,hfactfile,idisk1,iprint,id,nprocs,ierr,heade
  endif
 #endif
 
+ if (periodic) then
+    call set_boundary(xmin,xmax,ymin,ymax,zmin,zmax)
+ endif
+
 #ifdef ISOTHERMAL
  isothermal = .true.
 #else
@@ -674,7 +680,7 @@ end subroutine read_dump
 !+
 !-------------------------------------------------------------------
 subroutine read_smalldump(dumpfile,tfile,hfactfile,idisk1,iprint,id,nprocs,ierr,headeronly,dustydisc)
- use boundary,       only:xmin,xmax,ymin,ymax,zmin,zmax
+ use boundary,       only:set_boundary
  use dim,            only:maxp,gravity,maxalpha,mhd,use_dust,use_dustgrowth, &
                           h2chemistry,store_temperature,nsinkproperties
  use eos,            only:ieos,polyk,gamma,polyk2,qfacdisc,isink
@@ -687,7 +693,7 @@ subroutine read_smalldump(dumpfile,tfile,hfactfile,idisk1,iprint,id,nprocs,ierr,
                           ndustsmall,grainsize,graindens,Bextx,Bexty,Bextz, &
                           alphaind,poten,Bxyz,Bevol,dustfrac,deltav,        &
                           dustprop,tstop,VrelVf,temperature,abundance,      &
-                          ndusttypes,dustgasprop
+                          ndusttypes,dustgasprop,periodic
 #ifdef IND_TIMESTEPS
  use part,           only:dt_in
 #endif
@@ -709,6 +715,7 @@ subroutine read_smalldump(dumpfile,tfile,hfactfile,idisk1,iprint,id,nprocs,ierr,
 
  real(kind=4) :: dtind(npart)
 
+ real :: xmin,xmax,ymin,ymax,zmin,zmax
  character(len=200) :: fileident
  integer :: errors(5)
  logical :: smalldump,isothermal,ind_timesteps,const_av
@@ -799,6 +806,10 @@ subroutine read_smalldump(dumpfile,tfile,hfactfile,idisk1,iprint,id,nprocs,ierr,
  call allocate_memory(int(npart / nprocs))
 #endif
 
+ if (periodic) then
+    call set_boundary(xmin,xmax,ymin,ymax,zmin,zmax)
+ endif
+
 #ifdef ISOTHERMAL
  isothermal = .true.
 #else

src/main/writeheader.F90

@@ -155,9 +155,9 @@ subroutine write_header(icall,infile,evfile,logfile,dumpfile,ntot)
           write(iprint,"(2x,2(a,es14.6))") 'ymin = ',ymin,' ymax = ',ymax
           write(iprint,"(2x,2(a,es14.6))") 'zmin = ',zmin,' zmax = ',zmax
        else
-          write(iprint,"(2x,2(a,f10.6))")  'xmin = ',xmin,' xmax = ',xmax
-          write(iprint,"(2x,2(a,f10.6))")  'ymin = ',ymin,' ymax = ',ymax
-          write(iprint,"(2x,2(a,f10.6))")  'zmin = ',zmin,' zmax = ',zmax
+          write(iprint,"(2x,2(a,f10.5))")  'xmin = ',xmin,' xmax = ',xmax
+          write(iprint,"(2x,2(a,f10.5))")  'ymin = ',ymin,' ymax = ',ymax
+          write(iprint,"(2x,2(a,f10.5))")  'zmin = ',zmin,' zmax = ',zmax
        endif
     else
        write(iprint,"(a)") ' No boundaries set '

src/setup/phantomsetup.F90

@@ -87,12 +87,6 @@ program phantomsetup
  infile = trim(fileprefix)//'.in'
  inquire(file=trim(infile),exist=iexist)
 
- call set_default_options
-!
-!--if input file exists, read it
-!
- if (iexist) call read_infile(infile,logfile,evfile,dumpfile)
-
 !
 !--In general, setup routines do not know the number of particles until they
 !  are written. Need to allocate up to the hard limit. Legacy setup routines may
@@ -101,6 +95,12 @@ program phantomsetup
 !
  call allocate_memory(maxp_hard, part_only=.true.)
 
+ call set_default_options
+!
+!--if input file exists, read it
+!
+ if (iexist) call read_infile(infile,logfile,evfile,dumpfile)
+
 !
 !--reset logfile name
 !

src/setup/set_binary.f90

@@ -33,9 +33,12 @@ module setbinary
  use physcon, only:pi
  implicit none
  public :: set_binary,Rochelobe_estimate,L1_point,get_a_from_period
- public :: get_mean_angmom_vector
+ public :: get_mean_angmom_vector,get_eccentricity_vector
 
  private
+ interface get_eccentricity_vector
+  module procedure get_eccentricity_vector,get_eccentricity_vector_sinks
+ end interface get_eccentricity_vector
 
 contains
 
@@ -349,9 +352,9 @@ function get_eccentricity_vector(m1,m2,x1,x2,v1,v2)
  x0 = (m1*x1 + m2*x2)/(m1 + m2)
  v0 = (m1*v1 + m2*v2)/(m1 + m2)
 
- ! position and velocity vectors relative to centre of mass
- r = x2 - x0
- v = v2 - v0
+ ! position and velocity vectors relative to each other
+ r = x2 - x1
+ v = v2 - v1
 
  ! intermediate quantities
  dr = 1./sqrt(dot_product(r,r))
@@ -362,6 +365,25 @@ function get_eccentricity_vector(m1,m2,x1,x2,v1,v2)
 
 end function get_eccentricity_vector
 
+!----------------------------------------------------
+! interface to above assuming two sink particles
+!----------------------------------------------------
+function get_eccentricity_vector_sinks(xyzmh_ptmass,vxyz_ptmass,i1,i2)
+ real,    intent(in) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:)
+ integer, intent(in) :: i1, i2
+ real :: get_eccentricity_vector_sinks(3)
+
+ if (i1 > 0 .and. i2 > 0) then
+    get_eccentricity_vector_sinks = get_eccentricity_vector(&
+        xyzmh_ptmass(4,i1),xyzmh_ptmass(4,i2),&
+        xyzmh_ptmass(1:3,i1),xyzmh_ptmass(1:3,i2),&
+        vxyz_ptmass(1:3,i1),vxyz_ptmass(1:3,i2))
+ else
+    get_eccentricity_vector_sinks = 0.
+ endif
+
+end function get_eccentricity_vector_sinks
+
 !-------------------------------------------------------------
 ! Function to find mean angular momentum vector from a list
 ! of positions and velocities