Merge upstream updates

data/neutronstar/ns-rdensity.tab

@@ -1,5 +1,6 @@
- #   0.42442000000000002        1.6666700000000001        1.4249734000000001     
- #          914
+# r, density
+#   0.42442000000000002        1.6666700000000001        1.4249734000000001     
+#          914
    0.0000000000000000        1.4281200250341688     
    1.0952902519167577E-003   1.4281143125578766     
    2.1905805038335154E-003   1.4280971752432505     

scripts/buildbot.sh

@@ -77,6 +77,7 @@ listofcomponents='main setup analysis utils';
 # get list of targets, components and setups to check
 #
 allsetups=`grep 'ifeq ($(SETUP)' $phantomdir/build/Makefile_setups | grep -v skip | cut -d, -f 2 | cut -d')' -f 1`
+#allsetups='star'
 setuparr=($allsetups)
 batchsize=$(( ${#setuparr[@]} / $nbatch + 1 ))
 offset=$(( ($batch-1) * $batchsize ))
@@ -197,11 +198,11 @@ check_phantomsetup ()
    # run phantomsetup up to 3 times to successfully create/rewrite the .setup file
    #
    infile="${prefix}.in"
-   ./phantomsetup $prefix $flags < myinput.txt > /dev/null;
-   ./phantomsetup $prefix $flags < myinput.txt > /dev/null;
+   ./phantomsetup $prefix $flags < mycleanin.txt > /dev/null;
+   ./phantomsetup $prefix $flags < mycleanin.txt > /dev/null;
    if [ -e "$prefix.setup" ]; then
       print_result "creates .setup file" $pass;
-      #test_setupfile_options "$prefix" "$prefix.setup" $infile;
+      test_setupfile_options "$prefix.setup" "$flags" "$setup" $infile;
    else
       print_result "no .setup file" $warn;
    fi
@@ -255,22 +256,33 @@ check_phantomsetup ()
 test_setupfile_options()
 {
    myfail=0;
-   setup=$1;
-   setupfile=$2;
-   infile=$3;
+   #"$prefix.setup" "$flags" "$setup" $infile
+   setupfile=$1;
+   flags=$2;
+   setup=$3;
+   infile=$4;
    range=''
-   if [ "X$setup"=="Xstar" ]; then
-      param='iprofile'
-      range='1 2 3 4 5 6 7'
+   if [ "X$setup" == "Xstar" ]; then
+      param='iprofile1'
+      range='0 1 2 3 4 5 6 7'
    fi
    for x in $range; do
        valstring="$param = $x"
-       echo "checking $valstring"
+       echo "checking $valstring for SETUP=$setup"
+       rm $setupfile;
+       ./phantomsetup $setupfile $flags < mycleanin.txt > /dev/null;
        sed "s/$param.*=.*$/$valstring/" $setupfile > ${setupfile}.tmp
-       cp ${setupfile}.tmp $setupfile
+       mv ${setupfile}.tmp $setupfile
+       if [ "X$x" == "X6" ]; then
+          #sed "s/ieos.*=.*$/ieos = 9/" $setupfile > ${setupfile}.tmp
+          #mv ${setupfile}.tmp $setupfile
+          sed "s/dist_unit.*=.*$/dist_unit = km/" $setupfile > ${setupfile}.tmp
+          mv ${setupfile}.tmp $setupfile
+       fi
+       echo $setupfile
        rm $infile
-       ./phantomsetup $setupfile < /dev/null > /dev/null;
-       ./phantomsetup $setupfile < /dev/null;
+       ./phantomsetup $setupfile $flags < /dev/null > /dev/null;
+       ./phantomsetup $setupfile $flags < /dev/null > /dev/null;
 
        if [ -e $infile ]; then
           print_result "successful phantomsetup with $valstring" $pass;
@@ -424,7 +436,7 @@ for setup in $listofsetups; do
          echo $setup >> $faillog;
       fi
       if [ -e $errorlogold ]; then
-         diff --unchanged-line-format="" --old-line-format="" --new-line-format="%L" $errorlogold $errorlog | tail -20 > warnings.tmp
+         diff $errorlogold $errorlog | tail -20 > warnings.tmp
          if [ -s warnings.tmp ]; then
             newwarn=1;
          else

src/main/apr.f90

@@ -24,10 +24,12 @@ module apr
 !   mpiforce, part, physcon, ptmass, quitdump, random, relaxem,
 !   timestep_ind, vectorutils
 !
+ use dim, only:use_apr
  implicit none
 
  public :: init_apr,update_apr,read_options_apr,write_options_apr
  public :: create_or_update_apr_clump
+ public :: use_apr
 
  ! default values for runtime parameters
  integer, public :: apr_max_in = 3
@@ -56,7 +58,7 @@ module apr
 !+
 !-----------------------------------------------------------------------
 subroutine init_apr(apr_level,ierr)
- use dim, only:maxp_hard
+ use dim,  only:maxp_hard
  use part, only:npart,massoftype,aprmassoftype
  use apr_region, only:set_apr_centre,set_apr_regions
  integer,         intent(inout) :: ierr
@@ -610,22 +612,26 @@ subroutine write_options_apr(iunit)
  call write_inopt(apr_max_in,'apr_max','number of additional refinement levels (3 -> 2x resolution)',iunit)
  call write_inopt(ref_dir,'ref_dir','increase (1) or decrease (-1) resolution',iunit)
  call write_inopt(apr_type,'apr_type','1: static, 2: moving sink, 3: create clumps',iunit)
- select case (apr_type)
 
+ select case (apr_type)
  case(2)
     call write_inopt(track_part,'track_part','number of sink to track',iunit)
-
  case default
     call write_inopt(apr_centre(1),'apr_centre(1)','centre of region x position',iunit)
     call write_inopt(apr_centre(2),'apr_centre(2)','centre of region y position',iunit)
     call write_inopt(apr_centre(3),'apr_centre(3)','centre of region z position',iunit)
-
  end select
+
  call write_inopt(apr_rad,'apr_rad','radius of innermost region',iunit)
  call write_inopt(apr_drad,'apr_drad','size of step to next region',iunit)
 
 end subroutine write_options_apr
 
+!-----------------------------------------------------------------------
+!+
+!  Find the closest neighbour to a particle (needs replacing)
+!+
+!-----------------------------------------------------------------------
 subroutine closest_neigh(i,next_door,rmin)
  use part, only:xyzh,npart
  integer, intent(in)  :: i

src/main/eos.f90

@@ -47,8 +47,9 @@ module eos
 !   eos_stratified, infile_utils, io, mesa_microphysics, part, physcon,
 !   units
 !
- use part, only:ien_etotal,ien_entropy,ien_type
- use dim,  only:gr
+ use part,          only:ien_etotal,ien_entropy,ien_type
+ use dim,           only:gr
+ use eos_gasradrec, only:irecomb
  implicit none
  integer, parameter, public :: maxeos = 23
  real,               public :: polyk, polyk2, gamma
@@ -65,6 +66,8 @@ module eos
  public  :: init_eos,finish_eos,write_options_eos,read_options_eos
  public  :: write_headeropts_eos, read_headeropts_eos
 
+ public :: irecomb  ! propagated from eos_gasradrec
+
  private
 
  integer, public :: ieos          = 1
@@ -499,9 +502,8 @@ subroutine init_eos(eos_type,ierr)
     !--Check that if using ieos=6, then isink is set properly
     !
     if (isink==0) then
-       call error('eos','ieos=6, but isink is not set')
-       ierr = ierr_isink_not_set
-       return
+       call error('eos','ieos=6, but isink is not set, setting to 1')
+       isink = 1
     endif
 
  case(8)

src/main/evolve.F90

@@ -333,12 +333,9 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag)
        if (use_regnbody) then
           call group_identify(nptmass,n_group,n_ingroup,n_sing,xyzmh_ptmass,vxyz_ptmass,group_info,bin_info,nmatrix,&
                               new_ptmass=.true.,dtext=dtextforce)
-          call get_force(nptmass,npart,0,1,time,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass,vxyz_ptmass,&
-                 fxyz_ptmass,dsdt_ptmass,0.,0.,dummy,.false.,linklist_ptmass,bin_info,group_info=group_info)
-       else
-          call get_force(nptmass,npart,0,1,time,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass,vxyz_ptmass,&
-               fxyz_ptmass,dsdt_ptmass,0.,0.,dummy,.false.,linklist_ptmass,bin_info)
        endif
+       call get_force(nptmass,npart,0,1,time,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass,vxyz_ptmass,&
+                      fxyz_ptmass,dsdt_ptmass,0.,0.,dummy,.false.,linklist_ptmass,bin_info,group_info)
        if (ipart_createseeds /= 0) ipart_createseeds = 0 ! reset pointer to zero
        if (ipart_createstars /= 0) ipart_createstars = 0 ! reset pointer to zero
        dummy = 0

src/main/extern_densprofile.f90

@@ -16,7 +16,7 @@ module extern_densprofile
 !
 ! :Runtime parameters: None
 !
-! :Dependencies: datafiles, io, physcon, units
+! :Dependencies: datafiles, io, physcon, table_utils, units
 !
  implicit none
 
@@ -25,7 +25,7 @@ module extern_densprofile
  integer, public                         :: ntab
 
  ! *** Add option to .in file to specify density profile / mass enclosed filename? ***
- character(len=*), public, parameter     :: rhotabfile = 'density-profile.tab'
+ character(len=*), public, parameter     :: rhotabfile = 'relax1.profile'
  integer, public, parameter              :: nrhotab = 5000  ! maximum allowed size of r rho tabulated arrays
 
  public :: densityprofile_force, load_extern_densityprofile, read_rhotab, write_rhotab, calc_menc
@@ -41,6 +41,7 @@ module extern_densprofile
 !+
 !----------------------------------------------
 subroutine densityprofile_force(xi,yi,zi,fxi,fyi,fzi,phi)
+ use table_utils, only:yinterp
  real, intent(in)  :: xi, yi, zi
  real, intent(out) :: fxi, fyi, fzi, phi
 
@@ -58,7 +59,7 @@ subroutine densityprofile_force(xi,yi,zi,fxi,fyi,fzi,phi)
  if (ri2 >= r2surf) then
     f = fsurf * (r2surf / ri2)**1.5
  else
-    f = yinterp(ri2, r2tab, ftab, ntab)
+    f = yinterp(ftab(1:ntab),r2tab(1:ntab),ri2)
  endif
 
  fxi = f * xi
@@ -136,13 +137,15 @@ subroutine read_rhotab(filename, rsize, rtab, rhotab, nread, polyk, gamma, rhoc,
     return
  endif
 
- ! First line: # K gamma rhoc
+ ! first line, skip header
+ read(iunit, *,iostat=ierr)
+ ! second line: # K gamma rhoc
  read(iunit, *,iostat=ierr) hash,polyk, gamma, rhoc
  if (ierr /= 0) then
     call error('extern_densityprofile','Error reading first line of header from '//trim(filename))
     return
  endif
- ! Second line: # nentries  (number of r density entries in file)
+ ! third line: # nentries  (number of r density entries in file)
  read(iunit,*,iostat=ierr) hash,nread
  if (ierr /= 0) then
     call error('extern_densityprofile','Error reading second line of header from '//trim(filename))
@@ -180,11 +183,14 @@ subroutine write_rhotab(filename, rtab, rhotab, ntab, polyk, gamma, rhoc, ierr)
  ierr = 0
  open(newunit=iunit,file=filename,action='write',status='replace')
 
+ ! write header '# r,density'
+ write(iunit,"(a)") '# r, density'
+
  ! First line: # K gamma rhoc
- write(iunit,*) '# ', polyk, gamma, rhoc
+ write(iunit,"(a,3(g0,1x))") '# ', polyk, gamma, rhoc
 
  ! Second line: # nentries  (number of r density entries in file)
- write(iunit,*) '# ', ntab
+ write(iunit,"(a,i0)") '# ', ntab
 
  ! Loop over 'n' lines: r and density separated by space
  do i = 1,ntab
@@ -207,82 +213,17 @@ subroutine calc_menc(n, r, rho, menc_out, totmass)
  r2 = r(1:n)**2
  r2rho = r2(1:n) * rho(1:n)
 
- if (.false.) then
-    ! NB: Ensure that this mass calculation is correct if it is to be used (J. Wurster)
-    ! Use trapezoid, Simpson's and Simpson's 3/8 for first entries then Simpson's for remaining
-    ! (trapezoidal term has largest order error: avoid using it as part of sum for later terms)
-    menc(1) = 0.
-    menc(2) = (r(2)-r(1)) * (r2rho(1) + r2(2) * rho(2)) / 2.
-    menc(3) = (r(3)-r(1)) * (r2rho(1) + 4.*r2rho(2) + r2rho(3)) / 6.
-    menc(4) = (r(4)-r(1)) * (r2rho(1) + 3.*r2rho(2) + 3.*r2rho(3) + r2rho(4)) / 8.
-    totalmass = menc(1) + menc(2) + menc(3) + menc(4)
-    do i = 5, n
-       menc(i) = menc(i-2) + (r2rho(i-2) + 4.*r2rho(i-1) + r2rho(i)) * (r(i) - r(i-2)) / 6.
-       totalmass = totalmass + menc(i)
-    enddo
-    menc(:)   = 4.0 * pi * menc(:)
-    totalmass = 4.0 * pi * totalmass
- else
-    menc(1)   = 4.0/3.0*pi*r(1)**3 * rho(1)
-    do i = 2,n
-       menc(i)  = menc(i-1) + 4.0/3.0*pi*(r(i)**3 - r(i-1)**3) * rho(i)
-    enddo
-    totalmass = menc(n)
- endif
+ menc(1)   = 4.0/3.0*pi*r(1)**3 * rho(1)
+ do i = 2,n
+    menc(i)  = menc(i-1) + 4.0/3.0*pi*(r(i)**3 - r(i-1)**3) * rho(i)
+ enddo
+ totalmass = menc(n)
+
  if (present(menc_out)) menc_out = menc
  if (present(totmass))  totmass  = totalmass
 
 end subroutine calc_menc
 
-! Linear 1D interpolation
-real function yinterp(x, xtab, ytab, ntab)
- real,    intent(in) :: x
- real,    intent(in) :: ytab(:),xtab(:)
- integer, intent(in) :: ntab
-
- integer :: ibelow, iabove
- real    :: slope
-
- yinterp = 0.
- if (x <= xtab(1)) then
-    yinterp = ytab(1)
-    return
- elseif (x >= xtab(ntab)) then
-    yinterp = ytab(ntab)
-    return
- endif
-
- ibelow = indexbelow(x, xtab, ntab)
- iabove = ibelow + 1
-
- slope = (ytab(iabove) - ytab(ibelow)) / (xtab(iabove) - xtab(ibelow))
- yinterp = ytab(ibelow) + (x - xtab(ibelow)) * slope
-
-end function yinterp
-
-! Find index below value x in monotomic array xtab
-pure integer function indexbelow(x, xtab, ntab)
- real,    intent(in) :: x
- real,    intent(in) :: xtab(:)
- integer, intent(in) :: ntab
-
- integer :: ibelow, imid, iabove
-
- ibelow = 1
- iabove = ntab
- do while ((iabove - ibelow) > 1)
-    imid = (iabove + ibelow) / 2
-    if ((x > xtab(imid)) .eqv. (xtab(ntab) > xtab(1))) then
-       ibelow = imid
-    else
-       iabove = imid
-    endif
- enddo
-
- indexbelow = ibelow
-
-end function indexbelow
-
 !----------------------------------------------
 !+
 ! Wrapper to get the density profile (J.Wurster)

src/main/initial.F90

@@ -529,14 +529,10 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
     if (use_regnbody) then
        call init_subgroup
        call group_identify(nptmass,n_group,n_ingroup,n_sing,xyzmh_ptmass,vxyz_ptmass,group_info,bin_info,nmatrix)
-       call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,dtsinksink,&
-                             iexternalforce,time,merge_ij,merge_n,dsdt_ptmass,&
-                             group_info=group_info,bin_info=bin_info)
-
-    else
-       call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,dtsinksink,&
-                             iexternalforce,time,merge_ij,merge_n,dsdt_ptmass)
     endif
+    call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,dtsinksink,&
+                             iexternalforce,time,merge_ij,merge_n,dsdt_ptmass,&
+                             group_info,bin_info)
     dtsinksink = C_force*dtsinksink
     if (id==master) write(iprint,*) 'dt(sink-sink) = ',dtsinksink
     dtextforce = min(dtextforce,dtsinksink)
@@ -555,14 +551,9 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
           elseif (use_apr) then
              pmassi = aprmassoftype(igas,apr_level(i))
           endif
-          if (use_regnbody) then
-             call get_accel_sink_gas(nptmass,xyzh(1,i),xyzh(2,i),xyzh(3,i),xyzh(4,i),xyzmh_ptmass, &
-                                     fext(1,i),fext(2,i),fext(3,i),poti,pmassi,fxyz_ptmass,&
-                                     dsdt_ptmass,fonrmax,dtphi2,bin_info=bin_info)
-          else
-             call get_accel_sink_gas(nptmass,xyzh(1,i),xyzh(2,i),xyzh(3,i),xyzh(4,i),xyzmh_ptmass, &
-            fext(1,i),fext(2,i),fext(3,i),poti,pmassi,fxyz_ptmass,dsdt_ptmass,fonrmax,dtphi2)
-          endif
+          call get_accel_sink_gas(nptmass,xyzh(1,i),xyzh(2,i),xyzh(3,i),xyzh(4,i),xyzmh_ptmass, &
+                                  fext(1,i),fext(2,i),fext(3,i),poti,pmassi,fxyz_ptmass,&
+                                  dsdt_ptmass,fonrmax,dtphi2,bin_info)
           dtsinkgas = min(dtsinkgas,C_force*1./sqrt(fonrmax),C_force*sqrt(dtphi2))
        endif
     enddo