(setup) bug fixed in setup; use closepacked lattice; nptot no longer …

…optional in call to set_sphere to avoid seg fault

src/setup/set_sphere.f90

@@ -52,20 +52,20 @@ module spherical
 !  on equally spaced radial bins between rmin and rmax
 !+
 !-----------------------------------------------------------------------
-subroutine set_sphere(lattice,id,master,rmin,rmax,delta,hfact,np,xyzh, &
-                      rhofunc,rhotab,rtab,xyz_origin,nptot,dir,exactN,np_requested,mask)
+subroutine set_sphere(lattice,id,master,rmin,rmax,delta,hfact,np,xyzh,nptot, &
+                      rhofunc,rhotab,rtab,xyz_origin,dir,exactN,np_requested,mask)
  use stretchmap, only:set_density_profile,rho_func
  character(len=*), intent(in)    :: lattice
  integer,          intent(in)    :: id,master
  integer,          intent(inout) :: np
  real,             intent(in)    :: rmin,rmax,hfact
  real,             intent(out)   :: xyzh(:,:)
  real,             intent(inout) :: delta
+ integer(kind=8),  intent(inout) :: nptot
  procedure(rho_func), pointer, optional :: rhofunc
  real,             intent(in),    optional :: rhotab(:), rtab(:)
  integer,          intent(in),    optional :: dir
  integer,          intent(in),    optional :: np_requested
- integer(kind=8),  intent(inout), optional :: nptot
  real,             intent(in),    optional :: xyz_origin(3)
  logical,          intent(in),    optional :: exactN
  procedure(mask_prototype), optional :: mask
@@ -154,7 +154,7 @@ subroutine set_sphere_mc(id,master,rmin,rmax,hfact,np_requested,np,xyzh, &
  real,             intent(in)    :: rmin,rmax,hfact
  real,             intent(out)   :: xyzh(:,:)
  integer,          intent(out)   :: ierr
- integer(kind=8),  intent(inout), optional :: nptot
+ integer(kind=8),  intent(inout) :: nptot
  procedure(mask_prototype) :: mask
  integer :: i,npin,iseed,maxp
  real    :: vol_sphere,rr,phi,theta,mr,dir(3)
@@ -214,7 +214,7 @@ subroutine set_sphere_mc(id,master,rmin,rmax,hfact,np_requested,np,xyzh, &
     endif
  enddo
  ierr = 0
- if (present(nptot)) nptot = iparttot
+ nptot = iparttot
  if (id==master) write(*,"(1x,a,i10,a)") 'placed ',np-npin,' particles in random-but-symmetric sphere'
 
 end subroutine set_sphere_mc
@@ -394,6 +394,7 @@ subroutine set_unifdis_sphereN(lattice,id,master,xmin,xmax,ymin,ymax,zmin,zmax,p
  psep = (xmax-xmin)/nint((xmax-xmin)/psep)
 
 end subroutine set_unifdis_sphereN
+
 !-----------------------------------------------------------------------
 !+
 !  Wrapper to set an ellipse
@@ -431,5 +432,5 @@ subroutine set_ellipse(lattice,id,master,r_ellipsoid,delta,hfact,xyzh,np,nptot,n
  endif
 
 end subroutine set_ellipse
-!-----------------------------------------------------------------------
+
 end module spherical

src/setup/setup_solarsystem.f90

@@ -22,10 +22,9 @@ module setup
  implicit none
  public :: setpart
 
- real :: norbits
- integer :: dumpsperorbit,np_apophis
+ integer :: np_apophis
  logical :: asteroids
- character(len=20) :: epoch
+ character(len=20) :: epoch,tmax_in,dtmax_in
 
  private
 
@@ -39,7 +38,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
  use part,         only:nptmass,xyzmh_ptmass,vxyz_ptmass,idust,set_particle_type,&
                         grainsize,graindens,ndustlarge,ndusttypes,ndustsmall,ihacc,igas
  use setbinary,     only:set_binary
- use units,         only:set_units,umass,udist,unit_density,unit_velocity,utime
+ use units,         only:set_units,umass,udist,unit_density,unit_velocity,utime,in_code_units
  use physcon,       only:solarm,au,pi,km,solarr,ceresm,earthm,earthr,days
  use io,            only:master,fatal
  use timestep,      only:tmax,dtmax
@@ -49,6 +48,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
  use eos_tillotson, only:rho_0,A
  use spherical,     only:set_sphere
  use options,       only:ieos
+ use setup_params,  only:npart_total
  integer,           intent(in)    :: id
  integer,           intent(inout) :: npart
  integer,           intent(out)   :: npartoftype(:)
@@ -67,8 +67,8 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
 !
 ! default runtime parameters
 !
- norbits       = 1000.
- dumpsperorbit = 1
+ tmax_in = '1000 yr'
+ dtmax_in = '1 yr'
  asteroids = .true.
  np_apophis = 0
  !call date_and_time(values=values)
@@ -106,6 +106,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
 !--space available for injected gas particles
 !
  npart = 0
+ npart_total = 0
  npartoftype(:) = 0
  xyzh(:,:)  = 0.
  vxyzu(:,:) = 0.
@@ -115,8 +116,10 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
  mtot   = solarm/umass !  mass around which all bodies should orbit
 
  period = 2.*pi*sqrt(semia**3/mtot)
- tmax   = norbits*period
- dtmax  = period/dumpsperorbit
+ tmax   = in_code_units(tmax_in,ierr,unit_type='time')
+ if (ierr /= 0) call fatal('setup_solarsystem',' could not parse tmax')
+ dtmax  = in_code_units(dtmax_in,ierr,unit_type='time')
+ if (ierr /= 0) call fatal('setup_solarsystem',' could not parse dtmax')
 
  if (asteroids) then
     call set_minor_planets(npart,npartoftype,massoftype,xyzh,vxyzu,&
@@ -153,8 +156,11 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
     print "(3(a,1pg10.3),a)",' r_tidal is ',rtidal,' au,',rtidal*udist/km,' km, or ',rtidal*udist/earthr,' earth radii'
 
     if (np_apophis > 1) then
+       !
+       ! replace the sink particle with a ball of stuff
+       !
        dx = r_apophis/40.
-       call set_sphere('random',id,master,0.,r_apophis,dx,hfact,npart,xyzh,&
+       call set_sphere('closepacked',id,master,0.,r_apophis,dx,hfact,npart,xyzh,npart_total,&
                        xyz_origin=xyzmh_ptmass(1:3,nptmass),exactN=.true.,np_requested=np_apophis)
 
        do i=1,npart
@@ -163,10 +169,13 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
        massoftype(igas) = m_apophis / npart
        npartoftype(igas) = npart
        nptmass = nptmass - 1
+       !
+       ! print quantities from the equation of state to give an idea of the timestep
+       !
        if (ieos==23) then
-          spsoundmin = sqrt(A/rho_0)
-          print*,'     sound speed min = ',spsoundmin*unit_velocity,' cm/s'
-          print*,' sound crossing time = ',r_apophis/spsoundmin*utime/days,' days'
+          spsoundmin = sqrt(A/rho_0)/unit_velocity
+          print "(a,1pg11.4,a)",'     sound speed min = ',spsoundmin*unit_velocity/km,' km/s'
+          print "(a,1pg10.3,a)",' sound crossing time = ',(r_apophis/spsoundmin)*utime,' seconds'
        endif
     endif
  endif
@@ -193,8 +202,8 @@ subroutine write_setupfile(filename)
  open(unit=iunit,file=filename,status='replace',form='formatted')
 
  write(iunit,"(a)") '# input file for solar system setup routines'
- call write_inopt(norbits,'norbits','number of orbits at 1 au (i.e. end time in years)',iunit)
- call write_inopt(dumpsperorbit,'dumpsperorbit','number of dumps per orbit',iunit)
+ call write_inopt(tmax_in,'tmax_in','end time of simulation (e.g. 3 days)',iunit)
+ call write_inopt(dtmax_in,'dtmax_in','time between dumps (e.g. 1 hr)',iunit)
  call write_inopt(asteroids,'asteroids','add distant minor bodies as km-sized dust particles',iunit)
  call write_inopt(np_apophis,'np_apophis','number of particles used to represent apophis (0=none; 1=sink; n=gas)',iunit)
  call write_inopt(epoch,'epoch','epoch to query ephemeris, YYYY-MMM-DD HH:MM:SS.fff, blank = today',iunit)
@@ -219,8 +228,8 @@ subroutine read_setupfile(filename,ierr)
  nerr = 0
  ierr = 0
  call open_db_from_file(db,filename,iunit,ierr)
- call read_inopt(norbits,      'norbits',      db,min=0.,errcount=nerr)
- call read_inopt(dumpsperorbit,'dumpsperorbit',db,errcount=nerr)
+ call read_inopt(tmax_in, 'tmax_in',db,errcount=nerr)
+ call read_inopt(dtmax_in,'dtmax_in',db,errcount=nerr)
  call read_inopt(asteroids,'asteroids',db,errcount=nerr)
  call read_inopt(np_apophis,'np_apophis',db,errcount=nerr)
  call read_inopt(epoch,'epoch',db,errcount=nerr)