diff --git a/src/setup/set_sphere.f90 b/src/setup/set_sphere.f90 index e3358f9fd..8e3a72385 100644 --- a/src/setup/set_sphere.f90 +++ b/src/setup/set_sphere.f90 @@ -52,8 +52,8 @@ 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 @@ -61,11 +61,11 @@ subroutine set_sphere(lattice,id,master,rmin,rmax,delta,hfact,np,xyzh, & 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 diff --git a/src/setup/setup_solarsystem.f90 b/src/setup/setup_solarsystem.f90 index 5aafdb6f6..2aa99c84e 100644 --- a/src/setup/setup_solarsystem.f90 +++ b/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)