Merge pull request #416 from jameswurster/master

src/main/energies.F90

@@ -481,7 +481,7 @@ subroutine compute_energies(t)
                 do j = 9,21
                    n_ion = n_ion + data_out(j)
                 enddo
-                n_total = n_ion + data_out(4)
+                n_total = data_out(5)
                 if (n_total > 0.) then
                    n_total1 = 1.0/n_total
                 else
@@ -491,7 +491,7 @@ subroutine compute_energies(t)
                 call ev_data_update(ev_data_thread,iev_n(1),n_ion*n_total1)
                 call ev_data_update(ev_data_thread,iev_n(2),data_out( 8)*n_total1)
                 call ev_data_update(ev_data_thread,iev_n(3),data_out( 8))
-                call ev_data_update(ev_data_thread,iev_n(4),data_out( 4))
+                call ev_data_update(ev_data_thread,iev_n(4),n_total-n_ion)
                 call ev_data_update(ev_data_thread,iev_n(5),data_out(24))
                 call ev_data_update(ev_data_thread,iev_n(6),data_out(23))
                 call ev_data_update(ev_data_thread,iev_n(7),data_out(22))

src/setup/setup_sedov.f90

@@ -38,7 +38,7 @@ module setup
 !----------------------------------------------------------------
 subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact_out,time,fileprefix)
  use setup_params, only:rhozero
- use unifdis,      only:set_unifdis,latticetype,i_random,i_closepacked
+ use unifdis,      only:set_unifdis,latticetype,i_random,i_closepacked,get_xyzmin_xyzmax_exact
  use io,           only:iprint,master,fatal
  use boundary,     only:xmin,ymin,zmin,xmax,ymax,zmax,dxbound,dybound,dzbound
  use physcon,      only:pi
@@ -48,6 +48,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact_
  use part,         only:hfact,igas,periodic,set_particle_type
  use mpiutils,     only:bcast_mpi,reduceall_mpi
  use mpidomain,    only:i_belong
+ use prompting,    only:prompt
  use utils_shuffleparticles, only:shuffleparticles
  integer,           intent(in)    :: id
  integer,           intent(out)   :: npart
@@ -88,6 +89,13 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact_
        call fatal('setup','failed to read in all the data from .setup.  Aborting')
     endif
  elseif (id==master) then
+    print "(a,/)",trim(filename)//' not found: using interactive setup'
+    call prompt('Enter number of particles in x ',npartx,8,nint((maxp)**(1/3.)))
+    call prompt('Enter the type of particle lattice (1=cubic,2=closepacked,3=hcp,4=random)',ilattice,0,4)
+    shuffle_parts = .false.
+    if (ilattice==i_random) shuffle_parts = .true.
+    call prompt('Relax particles by shuffling?',shuffle_parts)
+
     call write_setupfile(filename)
     stop 'rerun phantomsetup after editing .setup file'
  else
@@ -108,6 +116,10 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact_
  gamma   = 5./3.
  gam1    = gamma - 1.
 
+ call get_xyzmin_xyzmax_exact(latticetype(ilattice),xmin,xmax,ymin,ymax,zmin,zmax,ierr,deltax)
+ dxbound = xmax-xmin
+ dybound = ymax-ymin
+ dzbound = zmax-zmin
  call set_unifdis(latticetype(ilattice),id,master,xmin,xmax,ymin,ymax,zmin,zmax,&
                   deltax,hfact,npart,xyzh,periodic,mask=i_belong)
 

src/setup/setup_unifdis.f90

@@ -58,7 +58,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
  use dim,          only:maxvxyzu,h2chemistry
  use setup_params, only:npart_total,ihavesetupB
  use io,           only:master
- use unifdis,      only:set_unifdis,latticetype
+ use unifdis,      only:set_unifdis,latticetype,get_xyzmin_xyzmax_exact
  use boundary,     only:xmin,ymin,zmin,xmax,ymax,zmax,dxbound,dybound,dzbound,set_boundary
  use part,         only:Bxyz,periodic,abundance,igas,iHI,dustfrac,ndustsmall,ndusttypes,grainsize,graindens
  use physcon,      only:pi,mass_proton_cgs,kboltz,years,pc,solarm,micron
@@ -199,6 +199,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
 
  call set_unifdis(latticetype(ilattice),id,master,xmin,xmax,ymin,ymax,zmin,zmax,&
                   deltax,hfact,npart,xyzh,periodic,nptot=npart_total,mask=i_belong)
+ call get_xyzmin_xyzmax_exact(latticetype(ilattice),xmin,xmax,ymin,ymax,zmin,zmax,ierr,deltax,npartx)
 
  npartoftype(:) = 0
  npartoftype(igas) = npart

src/utils/analysis_clumpfindWB23.F90

@@ -163,6 +163,8 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit)
     write(*,'(a,f7.2,a)')   '  The resolution of the ellipse-fitting routine to determine the clump size: ',dxgrid0_pc,' pc'
  endif
  write(*,'(a)')             '==============================================================='
+ write(*,'(a)')             ' WARNING! there is a bug fix on line labeled "Used in WB23"; the patch is untested'
+ write(*,'(a)')             '==============================================================='
  write(*,'(a)')             ' '
 
  !--calculate all densities