Skip to content

Commit

Permalink
Merge pull request #405 from danieljprice/set_star
Browse files Browse the repository at this point in the history
  • Loading branch information
@danieljprice
danieljprice committed Apr 21, 2023
2 parents 6e06a32 + dc3766b commit 6defda3
Show file tree
Hide file tree
Showing 13 changed files with 604 additions and 554 deletions.
26 changes: 11 additions & 15 deletions build/Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -637,7 +637,7 @@ OBJDUMP= $(OBJDUMP1:.F90=.o)
LIBSETUP=$(BINDIR)/libphantomsetup.a
SRCLIBSETUP=physcon.f90 geometry.f90 random.f90 utils_tables.f90 utils_vectors.f90 stretchmap.f90 \
utils_binary.f90 set_binary.f90 set_flyby.f90 \
set_hierarchical_utils.f90 set_hierarchical.f90 \
set_hierarchical_utils.f90 \
set_unifdis.f90 set_sphere.f90 set_shock.f90 \
set_dust.f90 libsetup.f90
OBJLIBSETUP=${SRCLIBSETUP:.f90=.o}
Expand All @@ -650,25 +650,21 @@ libsetup: $(OBJLIBSETUP)

#----------------------------------------------------
# these are the sources for phantom setup utility
# (just the list of extra files not included in libphantom)
#
.PHONY: phantomsetup
phantomsetup: setup

SRCSETUP= utils_summary.F90 utils_gravwave.f90 \
set_dust_options.f90 \
utils_indtimesteps.F90 partinject.F90 \
${SRCTURB} ${SRCNIMHD} ${SRCCHEM} \
ptmass.F90 energies.F90 density_profiles.f90 set_slab.f90 set_disc.F90 \
set_cubic_core.f90 set_fixedentropycore.f90 set_softened_core.f90 set_star.f90 relax_star.f90 \
set_vfield.f90 ptmass_radiation.F90 ${SRCINJECT} \
${SETUPFILE} checksetup.F90 \
set_Bfield.f90 readwrite_infile.f90 ${SRCKROME}
SRCSETUP= prompting.f90 set_dust_options.f90 \
density_profiles.f90 set_slab.f90 set_disc.F90 \
set_cubic_core.f90 set_fixedentropycore.f90 set_softened_core.f90 \
set_star_kepler.f90 set_star_mesa.f90 \
set_star.f90 relax_star.f90 set_hierarchical.f90 \
set_vfield.f90 set_Bfield.f90 \
${SETUPFILE}

OBJSETUP1= $(SRCSETUP:.f90=.o)
ifeq ($(KROME), krome)
OBJSETUP1 += $(KROME_OBJS)
endif
OBJSETUP= $(OBJDUMP) $(OBJSETUP1:.F90=.o) phantomsetup.o
OBJSETUP= $(OBJSETUP1:.F90=.o) phantomsetup.o

setup: checksystem checkparams libsetup libphantom $(OBJSETUP)
$(FC) $(FFLAGS) -o $(BINDIR)/phantomsetup $(OBJSETUP) $(LIBSETUP) $(LIBPHANTOM) $(LDFLAGS) $(LIBS)
Expand Down Expand Up @@ -800,7 +796,7 @@ interpolate3D_amr.o: adaptivemesh.o
ifndef MODDUMPBIN
MODDUMPBIN=phantommoddump
endif
OBJMOD1 = prompting.o set_Bfield.o density_profiles.o ${MODFILE:.f90=.o} phantom_moddump.o
OBJMOD1 = prompting.o set_Bfield.o density_profiles.o set_star_mesa.o ${MODFILE:.f90=.o} phantom_moddump.o
OBJMOD = ${OBJMOD1:.F90=.o}

phantom_moddump: checksystem checkparams libphantom libsetup $(OBJMOD)
Expand Down
73 changes: 39 additions & 34 deletions src/main/checksetup.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,15 +73,15 @@ subroutine check_setup(nerror,nwarn,restart)
endif

if (npart > maxp) then
print*,'Error in setup: npart (',npart,') > maxp (',maxp,')'
print*,'ERROR: npart (',npart,') > maxp (',maxp,')'
nerror = nerror + 1
endif
if (any(npartoftype < 0)) then
print*,'Error in setup: npartoftype -ve: ',npartoftype(:)
print*,'ERROR: npartoftype -ve: ',npartoftype(:)
nerror = nerror + 1
endif
if (sum(npartoftype) > maxp) then
print*,'Error in setup: sum(npartoftype) > maxp ',sum(npartoftype(:))
print*,'ERROR: sum(npartoftype) > maxp ',sum(npartoftype(:))
nerror = nerror + 1
endif
if (sum(npartoftype) /= npart) then
Expand All @@ -90,16 +90,16 @@ subroutine check_setup(nerror,nwarn,restart)
endif
#ifndef KROME
if (gamma <= 0.) then
print*,'WARNING! Error in setup: gamma not set (should be set > 0 even if not used)'
print*,'WARNING! gamma not set (should be set > 0 even if not used)'
nwarn = nwarn + 1
endif
#endif
if (hfact < 1. .or. hfact /= hfact) then
print*,'Error in setup: hfact = ',hfact,', should be >= 1'
print*,'ERROR: hfact = ',hfact,', should be >= 1'
nerror = nerror + 1
endif
if (polyk < 0. .or. polyk /= polyk) then
print*,'Error in setup: polyk = ',polyk,', should be >= 0'
print*,'ERROR: polyk = ',polyk,', should be >= 0'
nerror = nerror + 1
endif
#ifdef KROME
Expand All @@ -114,7 +114,7 @@ subroutine check_setup(nerror,nwarn,restart)
endif
#endif
if (npart < 0) then
print*,'Error in setup: npart = ',npart,', should be >= 0'
print*,'ERROR: npart = ',npart,', should be >= 0'
nerror = nerror + 1
elseif (npart==0 .and. nptmass==0) then
print*,'WARNING! setup: npart = 0 (and no sink particles either)'
Expand All @@ -130,12 +130,12 @@ subroutine check_setup(nerror,nwarn,restart)
!
if (all(iphase(1:npart)==0)) then
if (any(npartoftype(2:) > 0)) then
print*,'Error in setup: npartoftype > 0 for non-gas particles, but types have not been assigned'
print*,'ERROR: npartoftype > 0 for non-gas particles, but types have not been assigned'
nerror = nerror + 1
endif
iphase(1:npart) = isetphase(igas,iactive=.true.)
elseif (any(iphase(1:npart)==0)) then
print*,'Error in setup: types need to be assigned to all particles (or none)'
print*,'ERROR: types need to be assigned to all particles (or none)'
nerror = nerror + 1
endif
!
Expand Down Expand Up @@ -232,7 +232,7 @@ subroutine check_setup(nerror,nwarn,restart)
hmax = max(hi,hmax)
enddo
if (nbad > 0) then
print*,'Error in setup: negative, zero or ridiculous h on ',nbad,' of ',npart,' particles'
print*,'ERROR: negative, zero or ridiculous h on ',nbad,' of ',npart,' particles'
print*,' hmin = ',hmin,' hmax = ',hmax
nerror = nerror + 1
endif
Expand All @@ -249,12 +249,12 @@ subroutine check_setup(nerror,nwarn,restart)
endif
enddo
if (nbad > 0) then
print*,'Error in setup: negative thermal energy on ',nbad,' of ',npart,' particles'
print*,'ERROR: negative thermal energy on ',nbad,' of ',npart,' particles'
nerror = nerror + 1
endif
else
if (abs(gamma-1.) > tiny(gamma) .and. (ieos /= 2 .and. ieos /=9)) then
print*,'*** Error in setup: using isothermal EOS, but gamma = ',gamma
print*,'*** ERROR: using isothermal EOS, but gamma = ',gamma
gamma = 1.
print*,'*** Resetting gamma to 1, gamma = ',gamma
nwarn = nwarn + 1
Expand All @@ -269,7 +269,7 @@ subroutine check_setup(nerror,nwarn,restart)
nwarn = nwarn + 1
endif
if (npartoftype(itype) > 0 .and. .not.(in_range(massoftype(itype),0.))) then
print*,'Error in setup: massoftype = ',massoftype(itype),' for '//trim(labeltype(itype))// &
print*,'ERROR: massoftype = ',massoftype(itype),' for '//trim(labeltype(itype))// &
' particles (n'//trim(labeltype(itype))//' = ',npartoftype(itype),')'
nerror = nerror + 1
endif
Expand All @@ -279,7 +279,7 @@ subroutine check_setup(nerror,nwarn,restart)
!
call check_for_identical_positions(npart,xyzh,nbad)
if (nbad > 0) then
print*,'Error in setup: ',nbad,' of ',npart,' particles have identical or near-identical positions'
print*,'ERROR: ',nbad,' of ',npart,' particles have identical or near-identical positions'
nwarn = nwarn + 1
endif
!
Expand All @@ -296,7 +296,7 @@ subroutine check_setup(nerror,nwarn,restart)
endif
enddo
if (nbad > 0) then
print*,'Error in setup: ',nbad,' of ',npart,' particles setup OUTSIDE the periodic box'
print*,'ERROR: ',nbad,' of ',npart,' particles setup OUTSIDE the periodic box'
nerror = nerror + 1
endif
endif
Expand All @@ -317,15 +317,15 @@ subroutine check_setup(nerror,nwarn,restart)
if (accreted) nbad = nbad + 1
enddo
if (nbad > 0) then
print*,'Warning: ',nbad,' of ',npart,' particles setup within the accretion boundary'
print*,'WARNING: ',nbad,' of ',npart,' particles setup within the accretion boundary'
nwarn = nwarn + 1
endif
!--check if we are using a central accretor
hi = 0.
call accrete_particles(iexternalforce,0.,0.,0.,hi,massoftype(1),time,accreted)
!--if so, check for unresolved accretion radius
if (accreted .and. accradius1 < 0.5*hmin) then
print*,'Warning: accretion radius is unresolved by a factor of hmin/racc = ',hmin/accradius1
print*,'WARNING: accretion radius is unresolved by a factor of hmin/racc = ',hmin/accradius1
print*,'(this will cause the code to run needlessly slow)'
nwarn = nwarn + 1
endif
Expand All @@ -336,9 +336,9 @@ subroutine check_setup(nerror,nwarn,restart)
if (gravity .or. nptmass > 0) then
if (.not.G_is_unity()) then
if (gravity) then
print*,'Error in setup: self-gravity ON but G /= 1 in code units, got G=',get_G_code()
print*,'ERROR: self-gravity ON but G /= 1 in code units, got G=',get_G_code()
elseif (nptmass > 0) then
print*,'Error in setup: sink particles used but G /= 1 in code units, got G=',get_G_code()
print*,'ERROR: sink particles used but G /= 1 in code units, got G=',get_G_code()
endif
nerror = nerror + 1
endif
Expand All @@ -357,7 +357,7 @@ subroutine check_setup(nerror,nwarn,restart)
endif
if (mhd_nonideal) then
if (n_nden /= n_nden_phantom) then
print*,'Error in setup: n_nden in nicil.f90 needs to match n_nden_phantom in config.F90; n_nden = ',n_nden
print*,'ERROR: n_nden in nicil.f90 needs to match n_nden_phantom in config.F90; n_nden = ',n_nden
nerror = nerror + 1
endif
endif
Expand All @@ -367,7 +367,7 @@ subroutine check_setup(nerror,nwarn,restart)
!
if (npartoftype(idust) > 0) then
if (.not. use_dust) then
if (id==master) print*,'Error in setup: dust particles present but -DDUST is not set'
if (id==master) print*,'ERROR: dust particles present but -DDUST is not set'
nerror = nerror + 1
endif
if (use_dustfrac) then
Expand Down Expand Up @@ -422,22 +422,27 @@ subroutine check_setup(nerror,nwarn,restart)
!
call check_setup_ptmass(nerror,nwarn,hmin)
!
!--print centre of mass (must be done AFTER types have been checked)
!--check centre of mass
!
call get_centreofmass(xcom,vcom,npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass)
if (id==master) &
write(*,"(a,2(es10.3,', '),es10.3,a)") ' Centre of mass is at (x,y,z) = (',xcom,')'

if (.not.h2chemistry .and. maxvxyzu >= 4 .and. icooling == 3 .and. iexternalforce/=iext_corotate) then
if (dot_product(xcom,xcom) > 1.e-2) then
print*,'Error in setup: Gammie (2001) cooling (icooling=3) assumes Omega = 1./r^1.5'
print*,'ERROR: Gammie (2001) cooling (icooling=3) assumes Omega = 1./r^1.5'
print*,' but the centre of mass is not at the origin!'
nerror = nerror + 1
endif
endif

if (nerror==0 .and. nwarn==0) then
if (id==master) write(*,"(1x,a)") 'Particle setup OK'
!
!--print centre of mass (must be done AFTER types have been checked)
! ALSO, only print this if there are no warnings or errors to avoid obscuring warnings
!
if (id==master) then
write(*,"(a,2(es10.3,', '),es10.3,a)") ' Centre of mass is at (x,y,z) = (',xcom,')'
write(*,"(1x,a)") 'Particle setup OK'
endif
endif

end subroutine check_setup
Expand Down Expand Up @@ -486,17 +491,17 @@ subroutine check_setup_ptmass(nerror,nwarn,hmin)
real :: r,hsink

if (gr .and. nptmass > 0) then
print*,' Warning! Error in setup: nptmass = ',nptmass, ' should be = 0 for GR'
print*,' ERROR: nptmass = ',nptmass, ' should be = 0 for GR'
nwarn = nwarn + 1
return
endif

if (nptmass < 0) then
print*,' Error in setup: nptmass = ',nptmass, ' should be >= 0 '
print*,' ERROR: nptmass = ',nptmass, ' should be >= 0 '
nerror = nerror + 1
endif
if (nptmass > maxptmass) then
print*,' Error in setup: nptmass = ',nptmass,' exceeds ptmass array dimensions of ',maxptmass
print*,' ERROR: nptmass = ',nptmass,' exceeds ptmass array dimensions of ',maxptmass
nerror = nerror + 1
return
endif
Expand All @@ -512,10 +517,10 @@ subroutine check_setup_ptmass(nerror,nwarn,hmin)
dx = xyzmh_ptmass(1:3,j) - xyzmh_ptmass(1:3,i)
r = sqrt(dot_product(dx,dx))
if (r <= tiny(r)) then
print*,'Error in setup: sink ',j,' on top of sink ',i,' at ',xyzmh_ptmass(1:3,i)
print*,'ERROR: sink ',j,' on top of sink ',i,' at ',xyzmh_ptmass(1:3,i)
nerror = nerror + 1
elseif (r <= max(xyzmh_ptmass(ihacc,i),xyzmh_ptmass(ihacc,j))) then
print*,'Warning: sinks ',i,' and ',j,' within each others accretion radii: sep =',&
print*,'WARNING: sinks ',i,' and ',j,' within each others accretion radii: sep =',&
r,' h = ',xyzmh_ptmass(ihacc,i),xyzmh_ptmass(ihacc,j)
nwarn = nwarn + 1
endif
Expand All @@ -529,7 +534,7 @@ subroutine check_setup_ptmass(nerror,nwarn,hmin)
do i=1,nptmass
if (.not.in_range(xyzmh_ptmass(4,i))) then
nerror = nerror + 1
print*,' Error in setup: sink ',i,' mass = ',xyzmh_ptmass(4,i)
print*,' ERROR: sink ',i,' mass = ',xyzmh_ptmass(4,i)
elseif (xyzmh_ptmass(4,i) < 0.) then
print*,' Sink ',i,' has previously merged with another sink'
n = n + 1
Expand All @@ -544,11 +549,11 @@ subroutine check_setup_ptmass(nerror,nwarn,hmin)
hsink = max(xyzmh_ptmass(ihacc,i),xyzmh_ptmass(ihsoft,i))
if (hsink <= 0.) then
nerror = nerror + 1
print*,'Error in setup: sink ',i,' has accretion radius ',xyzmh_ptmass(ihacc,i),&
print*,'ERROR: sink ',i,' has accretion radius ',xyzmh_ptmass(ihacc,i),&
' and softening radius ',xyzmh_ptmass(ihsoft,i)
elseif (hsink <= 0.5*hmin .and. hmin > 0.) then
nwarn = nwarn + 1
print*,'Warning: sink ',i,' has unresolved accretion radius: hmin/racc = ',hmin/hsink
print*,'WARNING: sink ',i,' has unresolved accretion radius: hmin/racc = ',hmin/hsink
print*,' (this makes the code run pointlessly slow)'
endif
enddo
Expand Down
7 changes: 3 additions & 4 deletions src/main/io.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -218,7 +218,6 @@ subroutine formatreal(val,string,ierror,precision)
! endif
string = trim(adjustl(string))

return
end subroutine formatreal

!--------------------------------------------------------------------
Expand Down Expand Up @@ -406,12 +405,12 @@ subroutine warn(wherefrom,string,severity)
if (buffer_warnings) then
call buffer_warning(trim(wherefrom),trim(string),ncount)
if (ncount < maxcount) then
write(iprint,"(' WARNING! ',a,': ',a)") trim(wherefrom),trim(string)
write(iprint,"(/' WARNING! ',a,': ',a,/)") trim(wherefrom),trim(string)
elseif (ncount==maxcount) then
write(iprint,"(' (buffering remaining warnings... ',a,') ')") trim(string)
endif
else
write(iprint,"(' WARNING! ',a,': ',a)") trim(wherefrom),trim(string)
write(iprint,"(/' WARNING! ',a,': ',a,/)") trim(wherefrom),trim(string)
endif
case(3)
ncount = 0_8
Expand All @@ -434,7 +433,6 @@ subroutine warn(wherefrom,string,severity)
write(iprint,"(/' WARNING(unknown severity)! ',a,': ',a)") trim(wherefrom),trim(string)
end select

return
end subroutine warn

!--------------------------------------------------------------------
Expand Down Expand Up @@ -565,6 +563,7 @@ subroutine die
#endif

call exit(1)

end subroutine die

end module io
Loading

0 comments on commit 6defda3

Please sign in to comment.