Relax bug in gr

.github/workflows/krome.yml

@@ -33,6 +33,10 @@ jobs:
         with:
           compiler: ${{ matrix.toolchain.compiler }}
 
+      - name: "Install LAPACK dependencies (Ubuntu)"
+        if: matrix.os == 'ubuntu-latest'
+        run: sudo apt-get update && sudo apt-get install -y liblapack-dev
+
       - name: "Clone phantom"
         uses: actions/checkout@v4
 

build/Makefile

@@ -1179,20 +1179,20 @@ evol_dustywaves: checksys $(OBJDUSTEVOL)
 # these are the sources for the ev2mdot utility
 #
 .PHONY: ev2mdot
-ev2mdot: checksys utils_filenames.o utils_infiles.o utils_evfiles.o ev2mdot.o
+ev2mdot: checksys utils_filenames.o utils_infiles.o utils_evfiles.o ev2dot_utils.o ev2mdot.o
 	@echo ""
 	@echo "ev2mdot: Accretion rates R us."
 	@echo ""
-	$(FC) $(FFLAGS) -o $(BINDIR)/$@ ev2mdot.o utils_filenames.o utils_evfiles.o utils_infiles.o
+	$(FC) $(FFLAGS) -o $(BINDIR)/$@ ev2mdot.o utils_filenames.o utils_infiles.o utils_evfiles.o ev2dot_utils.o
 
 cleanev2mdot:
 	rm -f $(BINDIR)/ev2mdot
 #----------------------------------------------------
 # these are the sources for the ev2kdot utility
 #
 .PHONY: ev2kdot
-ev2kdot: checksys utils_filenames.o utils_infiles.o utils_evfiles.o ev2kdot.o
-	$(FC) $(FFLAGS) -o $(BINDIR)/$@ ev2kdot.o utils_filenames.o utils_evfiles.o utils_infiles.o
+ev2kdot: checksys utils_filenames.o utils_infiles.o utils_evfiles.o ev2dot_utils.o ev2kdot.o
+	$(FC) $(FFLAGS) -o $(BINDIR)/$@ ev2kdot.o utils_filenames.o utils_infiles.o utils_evfiles.o ev2dot_utils.o
 
 cleanev2kdot:
 	rm -f $(BINDIR)/ev2kdot

scripts/writemake.sh

@@ -12,7 +12,7 @@ echo '#'
 echo 'PHANTOMDIR='${0/scripts\/writemake.sh/};
 echo 'SPLASHDIR='$splashdir
 echo 'EDITOR=vi'
-makeflags='RUNDIR=${PWD}';
+makeflags='RUNDIR="${PWD}"';
 if [ $# -ge 1 ]; then
    echo 'ifndef SETUP';
    echo 'SETUP='$1;

src/main/eos.f90

@@ -1498,7 +1498,7 @@ end subroutine write_headeropts_eos
 !-----------------------------------------------------------------------
 subroutine read_headeropts_eos(ieos,hdr,ierr)
  use dump_utils,        only:dump_h, extract
- use io,                only:iprint,id,master
+ use io,                only:iprint,id,master,iverbose
  use dim,               only:use_krome,maxvxyzu
  integer,      intent(in)  :: ieos
  type(dump_h), intent(in)  :: hdr
@@ -1512,12 +1512,12 @@ subroutine read_headeropts_eos(ieos,hdr,ierr)
  if (id==master) then
     if (maxvxyzu >= 4) then
        if (use_krome) then
-          write(iprint,*) 'KROME eos: initial gamma = 1.666667'
+         if (iverbose >= 0) write(iprint,*) 'KROME eos: initial gamma = 1.666667'
        else
-          write(iprint,*) 'adiabatic eos: gamma = ',gamma
+          if (iverbose >= 0) write(iprint,*) 'adiabatic eos: gamma = ',gamma
        endif
     else
-       write(iprint,*) 'setting isothermal sound speed^2 (polyk) = ',polyk,' gamma = ',gamma
+      if (iverbose >= 0) write(iprint,*) 'setting isothermal sound speed^2 (polyk) = ',polyk,' gamma = ',gamma
        if (polyk <= tiny(polyk)) write(iprint,*) 'WARNING! sound speed zero in dump!, polyk = ',polyk
     endif
  endif
@@ -1536,7 +1536,7 @@ subroutine read_headeropts_eos(ieos,hdr,ierr)
        if (id==master) write(iprint,*) 'ERROR: qfacdisc <= 0'
        ierr = 2
     else
-       if (id==master) write(iprint,*) 'qfacdisc = ',qfacdisc
+       if (id==master .and. iverbose >= 0) write(iprint,*) 'qfacdisc = ',qfacdisc
     endif
  endif
 
@@ -1549,7 +1549,7 @@ subroutine read_headeropts_eos(ieos,hdr,ierr)
        if (id==master) write(iprint,*) 'ERROR: qfacdisc2 == 0'
        ierr = 2
     else
-       if (id==master) write(iprint,*) 'qfacdisc2 = ',qfacdisc2
+       if (id==master .and. iverbose >= 0) write(iprint,*) 'qfacdisc2 = ',qfacdisc2
     endif
  endif
 

src/main/readwrite_dumps_common.f90

@@ -120,7 +120,7 @@ subroutine unfill_header(hdr,phantomdump,got_tags,nparttot, &
                          nblocks,npart,npartoftype, &
                          tfile,hfactfile,alphafile,iprint,id,nprocs,ierr)
  use dim,        only:maxdustlarge,use_dust
- use io,         only:master ! check this
+ use io,         only:master,iverbose ! check this
  use eos,        only:isink
  use part,       only:maxtypes,igas,idust,ndustsmall,ndustlarge,ndusttypes,&
                       npartoftypetot
@@ -203,7 +203,7 @@ subroutine unfill_header(hdr,phantomdump,got_tags,nparttot, &
  if (nblocks==1) then
     npart = int(nparttoti)
     nparttot = npart
-    if (id==master) write (iprint,*) 'npart = ',npart
+    if (id==master .and. iverbose >= 0) write (iprint,*) 'npart = ',npart
  endif
  if (got_tags) then
     call extract('ntypes',ntypesinfile8,hdr,ierr1)
@@ -218,7 +218,7 @@ subroutine unfill_header(hdr,phantomdump,got_tags,nparttot, &
  if (nblocks > 1) then
     call extract('npartoftype',npartoftype(1:ntypesinfile),hdr,ierr1)
  endif
- if (id==master) write(*,*) 'npart(total) = ',nparttot
+ if (id==master .and. iverbose >= 0) write(*,*) 'npart(total) = ',nparttot
 !
 !--number of dust species
 !
@@ -250,7 +250,7 @@ subroutine unfill_header(hdr,phantomdump,got_tags,nparttot, &
                      tfile,hfactfile,alphafile,iprint,ierr)
  if (ierr /= 0) return
 
- if (id==master) write(iprint,*) 'time = ',tfile
+ if (id==master .and. iverbose >= 0) write(iprint,*) 'time = ',tfile
 
 end subroutine unfill_header
 
@@ -363,6 +363,7 @@ subroutine fill_header(sphNGdump,t,nparttot,npartoftypetot,nblocks,nptmass,hdr,i
        call add_to_rheader(rho_bkg_ini,'rho_bkg_ini',hdr,ierr)
     endif
     call add_to_rheader(get_conserv,'get_conserv',hdr,ierr)
+    call add_to_rheader(mtot_in,'mtot_in',hdr,ierr)
     call add_to_rheader(etot_in,'etot_in',hdr,ierr)
     call add_to_rheader(angtot_in,'angtot_in',hdr,ierr)
     call add_to_rheader(totmom_in,'totmom_in',hdr,ierr)
@@ -371,9 +372,6 @@ subroutine fill_header(sphNGdump,t,nparttot,npartoftypetot,nblocks,nptmass,hdr,i
        call add_to_rheader(grainsize(1:ndusttypes),'grainsize',hdr,ierr)
        call add_to_rheader(graindens(1:ndusttypes),'graindens',hdr,ierr)
     endif
-    if (use_apr) then
-       call add_to_rheader(mtot_in,'mtot_in',hdr,ierr)
-    endif
  endif
 
  ! real*8

src/main/readwrite_dumps_fortran.f90

@@ -524,7 +524,7 @@ subroutine read_dump_fortran(dumpfile,tfile,hfactfile,idisk1,iprint,id,nprocs,ie
  type(dump_h)          :: hdr
  integer               :: i,ierrh
 
- if (id==master) write(iprint,"(/,1x,a,i3)") '>>> reading setup from file: '//trim(dumpfile)//' on unit ',idisk1
+ if (id==master .and. iverbose >= 0) write(iprint,"(/,1x,a,i3)") '>>> reading setup from file: '//trim(dumpfile)//' on unit ',idisk1
  opened_full_dump_fortran = .true.
  dt_read_in_fortran       = .false.
  !
@@ -558,7 +558,7 @@ subroutine read_dump_fortran(dumpfile,tfile,hfactfile,idisk1,iprint,id,nprocs,ie
     close(idisk1)
     return
  endif
- if (id==master) write(iprint,*) trim(fileidentr)
+ if (id==master .and. iverbose >= 0) write(iprint,*) trim(fileidentr)
 
  ! extract file type from the fileid string
  call get_options_from_fileid(fileidentr,tagged,phantomdump,smalldump,use_dustfrac,ierr)
@@ -672,7 +672,7 @@ subroutine read_dump_fortran(dumpfile,tfile,hfactfile,idisk1,iprint,id,nprocs,ie
     elseif (npartread > 0) then
        string = ''
        if (nprocs > 1) write(string,'(a,i5)') 'thread',iblock
-       write(*,"(2(a,i10),a,i5,a,i10,'-',i10)") trim(string)//' reading particles ',noffset+1,&
+       if (iverbose >= 0) write(*,"(2(a,i10),a,i5,a,i10,'-',i10)") trim(string)//' reading particles ',noffset+1,&
            ':',noffset+npartread,', from block ',iblock,' lims=',i1,i2
     else
        write(*,"(a,i10,a)") ' WARNING! block contains no SPH particles, reading ',nptmass,' point mass particles only'
@@ -722,9 +722,9 @@ subroutine read_dump_fortran(dumpfile,tfile,hfactfile,idisk1,iprint,id,nprocs,ie
  endif
 
  if (narraylengths >= 4) then
-    if (id==master) write(iprint,"(a,/)") ' <<< finished reading (MHD) file '
+    if (id==master .and. iverbose >= 0) write(iprint,"(a,/)") ' <<< finished reading (MHD) file '
  else
-    if (id==master) write(iprint,"(a,/)") ' <<< finished reading (hydro) file '
+    if (id==master .and. iverbose >= 0) write(iprint,"(a,/)") ' <<< finished reading (hydro) file '
  endif
  close(idisk1)
  return
@@ -802,7 +802,7 @@ subroutine read_smalldump_fortran(dumpfile,tfile,hfactfile,idisk1,iprint,id,npro
     ierr = 1
     return
  endif
- if (id==master) write(iprint,*) trim(fileidentr)
+ if (id==master .and. iverbose >= 0) write(iprint,*) trim(fileidentr)
 
  ! extract file type from the fileid string
  call get_options_from_fileid(fileidentr,tagged,phantomdump,smalldump,use_dustfrac,ierr)

src/main/utils_sort.f90

@@ -657,10 +657,13 @@ subroutine sort_by_radius(n,xyzh,iorder,x0)
  real, intent(in), optional :: x0(3)
 
  ! optional argument x0=[1,1,1] to set the origin
- if (present(x0)) call set_r2func_origin(x0(1),x0(2),x0(3))
-
- ! sort by r^2 using the r2func function
- call indexxfunc(n,r2func,xyzh,iorder)
+ if (present(x0)) then
+    call set_r2func_origin(x0(1),x0(2),x0(3))
+    call indexxfunc(n,r2func_origin,xyzh,iorder)
+ else
+    ! sort by r^2 using the r2func function
+    call indexxfunc(n,r2func,xyzh,iorder)
+ endif
 
 end subroutine sort_by_radius
 

src/setup/readwrite_kepler.f90

@@ -199,10 +199,11 @@ end subroutine read_kepler_file
 !
 !+
 !-----------------------------------------------------------------------
-subroutine write_kepler_comp(composition,comp_label,columns_compo,r,&
+subroutine write_kepler_comp(filename,composition,comp_label,columns_compo,r,&
                              xyzh,npart,npts,composition_exists,npin)
 
  use table_utils, only:yinterp
+ character(len=*), intent(in)               :: filename
  integer, intent(in)                        :: columns_compo,npart,npts
  real,    intent(in)                        :: xyzh(:,:)
  real, allocatable,intent(in)               :: r(:)
@@ -212,7 +213,7 @@ subroutine write_kepler_comp(composition,comp_label,columns_compo,r,&
  integer, intent(in), optional              :: npin
  real , allocatable                         :: compositioni(:,:)
  real, allocatable                          :: comp(:)
- integer                                    :: i,j,iu,i1
+ integer                                    :: i,j,iu,i1,ierr
  real                                       :: ri
 
  composition_exists = .false.
@@ -226,9 +227,13 @@ subroutine write_kepler_comp(composition,comp_label,columns_compo,r,&
  endif
 
  if (composition_exists) then
-    print*, 'Writing the stellar composition for each particle into ','kepler.comp'
+    print*, 'Writing the stellar composition for each particle into '//trim(filename)
 
-    open(newunit=iu,file='kepler.comp')
+    open(newunit=iu,file=trim(filename),status='replace',form='formatted',iostat=ierr)
+    if (ierr /= 0) then
+       print*,' ERROR writing to '//trim(filename)
+       return
+    endif
     write(iu,"('#',50(1x,'[',1x,a7,']',2x))") &
             comp_label
     !Now setting the composition of star if the case used was ikepler
@@ -246,7 +251,6 @@ subroutine write_kepler_comp(composition,comp_label,columns_compo,r,&
             (compositioni(j,1),j=1,columns_compo)
     enddo
     close(iu)
-    print*, '>>>>>> done'
  endif
 
 end subroutine write_kepler_comp