Merge pull request #548 from Yrisch/fixes_and_patches

Fix the corruption issue while creating sink dynamically

src/main/evolve.F90

@@ -276,7 +276,7 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag)
        ! creation of new sink particles
        !
        call ptmass_create(nptmass,npart,ipart_rhomax,xyzh,vxyzu,fxyzu,fext,divcurlv,&
-                          poten,massoftype,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dptmass,time)
+                          poten,massoftype,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,fxyz_ptmass_sinksink,dptmass,time)
     endif
     !
     ! Strang splitting: implicit update for half step

src/main/initial.F90

@@ -127,7 +127,7 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
  use part,             only:npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,Bevol,dBevol,tau, tau_lucy, &
                             npartoftype,maxtypes,ndusttypes,alphaind,ntot,ndim,update_npartoftypetot,&
                             maxphase,iphase,isetphase,iamtype,igas,idust,imu,igamma,massoftype, &
-                            nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,&
+                            nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,fxyz_ptmass_sinksink,&
                             epot_sinksink,get_ntypes,isdead_or_accreted,dustfrac,ddustevol,&
                             nden_nimhd,dustevol,rhoh,gradh, &
                             Bevol,Bxyz,dustprop,filfac,ddustprop,ndustsmall,iboundary,eos_vars,dvdx, &
@@ -546,6 +546,9 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
     if (r_merge_uncond < 2.0*h_acc) then
        write(iprint,*) ' WARNING! Sink creation is on, but but merging is off!  Suggest setting r_merge_uncond >= 2.0*h_acc'
     endif
+    dsdt_ptmass = 0. ! could introduce NaN in ptmass spins if not initialised (no get_accel done before creating sink)
+    fxyz_ptmass = 0.
+    fxyz_ptmass_sinksink = 0.
  endif
  if (abs(time) <= tiny(0.)) then
     !initialize nucleation array at the start of the run only

src/main/ptmass.F90

@@ -752,7 +752,7 @@ subroutine ptmass_accrete(is,nptmass,xi,yi,zi,hi,vxi,vyi,vzi,fxi,fyi,fzi, &
                           dptmass,time,facc,nbinmax,ibin_wakei,nfaili)
 
 !$ use omputils, only:ipart_omp_lock
- use part,       only: ihacc
+ use part,       only: ihacc,ndptmass
  use kernel,     only: radkern2
  use io,         only: iprint,iverbose,fatal
  use io_summary, only: iosum_ptmass,maxisink,print_acc
@@ -762,7 +762,7 @@ subroutine ptmass_accrete(is,nptmass,xi,yi,zi,hi,vxi,vyi,vzi,fxi,fyi,fzi, &
  real,              intent(in)    :: xyzmh_ptmass(nsinkproperties,nptmass)
  real,              intent(in)    :: vxyz_ptmass(3,nptmass)
  logical,           intent(out)   :: accreted
- real,              intent(inout) :: dptmass(:,:)
+ real,              intent(inout) :: dptmass(ndptmass,nptmass)
  integer(kind=1),   intent(in)    :: nbinmax
  integer(kind=1),   intent(inout) :: ibin_wakei
  integer, optional, intent(out)   :: nfaili
@@ -936,11 +936,13 @@ end subroutine ptmass_accrete
 !+
 !-----------------------------------------------------------------------
 subroutine update_ptmass(dptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,nptmass)
- real,    intent(in)    :: dptmass(:,:)
+ use part ,only:ndptmass
+ integer, intent(in)    :: nptmass
+ real,    intent(in)    :: dptmass(ndptmass,nptmass)
  real,    intent(inout) :: xyzmh_ptmass(:,:)
  real,    intent(inout) :: vxyz_ptmass(:,:)
  real,    intent(inout) :: fxyz_ptmass(:,:)
- integer, intent(in)    :: nptmass
+
 
  real                   :: newptmass(nptmass),newptmass1(nptmass)
 
@@ -996,9 +998,9 @@ end subroutine update_ptmass
 !+
 !-------------------------------------------------------------------------
 subroutine ptmass_create(nptmass,npart,itest,xyzh,vxyzu,fxyzu,fext,divcurlv,poten,&
-                         massoftype,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dptmass,time)
+                         massoftype,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,fxyz_ptmass_sinksink,dptmass,time)
  use part,   only:ihacc,ihsoft,igas,iamtype,get_partinfo,iphase,iactive,maxphase,rhoh, &
-                  ispinx,ispiny,ispinz,fxyz_ptmass_sinksink,eos_vars,igasP,igamma,ndptmass
+                  ispinx,ispiny,ispinz,eos_vars,igasP,igamma,ndptmass
  use dim,    only:maxp,maxneigh,maxvxyzu,maxptmass,ind_timesteps
  use kdtree, only:getneigh
  use kernel, only:kernel_softening,radkern
@@ -1022,8 +1024,8 @@ subroutine ptmass_create(nptmass,npart,itest,xyzh,vxyzu,fxyzu,fext,divcurlv,pote
  real,            intent(inout) :: xyzh(:,:)
  real,            intent(in)    :: vxyzu(:,:),fxyzu(:,:),fext(:,:),massoftype(:)
  real(4),         intent(in)    :: divcurlv(:,:),poten(:)
- real,            intent(inout) :: xyzmh_ptmass(:,:)
- real,            intent(inout) :: vxyz_ptmass(:,:),fxyz_ptmass(:,:),dptmass(ndptmass,nptmass+1)
+ real,            intent(inout) :: xyzmh_ptmass(:,:),dptmass(ndptmass,maxptmass)
+ real,            intent(inout) :: vxyz_ptmass(:,:),fxyz_ptmass(4,maxptmass),fxyz_ptmass_sinksink(4,maxptmass)
  real,            intent(in)    :: time
  integer(kind=1)    :: iphasei,ibin_wakei,ibin_itest
  integer            :: nneigh
@@ -1533,8 +1535,8 @@ subroutine ptmass_create(nptmass,npart,itest,xyzh,vxyzu,fxyzu,fext,divcurlv,pote
     nacc = int(reduceall_mpi('+', nacc))
 
     ! update ptmass position, spin, velocity, acceleration, and mass
-    fxyz_ptmass(:,nptmass) = 0.0
-    fxyz_ptmass_sinksink(:,nptmass) = 0.0
+    fxyz_ptmass(1:4,n) = 0.0
+    fxyz_ptmass_sinksink(1:4,n) = 0.0
     call update_ptmass(dptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,nptmass)
 
     if (id==id_rhomax) then

src/main/substepping.F90

@@ -430,7 +430,7 @@ subroutine substep(npart,ntypes,nptmass,dtsph,dtextforce,time,xyzh,vxyzu,fext, &
                    n_group,n_ingroup,n_sing)
  use io,             only:iverbose,id,master,iprint,fatal
  use options,        only:iexternalforce
- use part,           only:fxyz_ptmass_sinksink
+ use part,           only:fxyz_ptmass_sinksink,ndptmass
  use io_summary,     only:summary_variable,iosumextr,iosumextt
  use externalforces, only:is_velocity_dependent
  use ptmass,         only:use_fourthorder,use_regnbody,ck,dk
@@ -442,7 +442,7 @@ subroutine substep(npart,ntypes,nptmass,dtsph,dtextforce,time,xyzh,vxyzu,fext, &
  real,            intent(inout) :: dtextforce
  real,            intent(inout) :: xyzh(:,:),vxyzu(:,:),fext(:,:)
  real,            intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:),fxyz_ptmass(:,:),dsdt_ptmass(:,:)
- real,            intent(inout) :: dptmass(:,:),fsink_old(:,:),gtgrad(:,:)
+ real,            intent(inout) :: dptmass(ndptmass,nptmass),fsink_old(:,:),gtgrad(:,:)
  integer(kind=1), intent(in)    :: nbinmax
  integer(kind=1), intent(inout) :: ibin_wake(:),nmatrix(nptmass,nptmass)
  logical :: extf_vdep_flag,done,last_step,accreted

src/main/utils_dumpfiles.f90

@@ -1358,56 +1358,56 @@ function allocate_header(nint,nint1,nint2,nint4,nint8,nreal,nreal4,nreal8,err) r
  integer, intent(in),  optional :: nint,nint1,nint2,nint4,nint8,nreal,nreal4,nreal8
  integer, intent(out), optional :: err
  type(dump_h) :: hdr
- integer      :: size(ndatatypes)
+ integer      :: size_type(ndatatypes)
  integer      :: ierrs(ndatatypes)
  integer      :: ierr
 
  ! make sure header is deallocated first
  call free_header(hdr,ierr)
 
- size(:) = maxphead
- if (present(nint))   size(i_int)  = nint
- if (present(nint1))  size(i_int1) = nint1
- if (present(nint2))  size(i_int2) = nint2
- if (present(nint4))  size(i_int4) = nint4
- if (present(nint8))  size(i_int8) = nint8
- if (present(nreal))  size(i_real) = nreal
- if (present(nreal4)) size(i_real4) = nreal4
- if (present(nreal8)) size(i_real8) = nreal8
+ size_type(:) = maxphead
+ if (present(nint))   size_type(i_int)  = nint
+ if (present(nint1))  size_type(i_int1) = nint1
+ if (present(nint2))  size_type(i_int2) = nint2
+ if (present(nint4))  size_type(i_int4) = nint4
+ if (present(nint8))  size_type(i_int8) = nint8
+ if (present(nreal))  size_type(i_real) = nreal
+ if (present(nreal4)) size_type(i_real4) = nreal4
+ if (present(nreal8)) size_type(i_real8) = nreal8
 
  if (present(err)) err = 0
  ierrs(:) = 0
  hdr%nums(:) = 0
- if (size(i_int) > 0)  then
-    allocate(hdr%inttags(size(i_int)),hdr%intvals(size(i_int)),stat=ierrs(1))
+ if (size_type(i_int) > 0)  then
+    allocate(hdr%inttags(size_type(i_int)),hdr%intvals(size_type(i_int)),stat=ierrs(1))
     if (ierrs(1)==0) hdr%inttags(:) = ''
  endif
- if (size(i_int1) > 0) then
-    allocate(hdr%int1tags(size(i_int1)),hdr%int1vals(size(i_int1)),stat=ierrs(2))
+ if (size_type(i_int1) > 0) then
+    allocate(hdr%int1tags(size_type(i_int1)),hdr%int1vals(size_type(i_int1)),stat=ierrs(2))
     if (ierrs(2)==0) hdr%int1tags(:) = ''
  endif
- if (size(i_int2) > 0) then
-    allocate(hdr%int2tags(size(i_int2)),hdr%int2vals(size(i_int2)),stat=ierrs(3))
+ if (size_type(i_int2) > 0) then
+    allocate(hdr%int2tags(size_type(i_int2)),hdr%int2vals(size_type(i_int2)),stat=ierrs(3))
     if (ierrs(3)==0) hdr%int2tags(:) = ''
  endif
- if (size(i_int4) > 0) then
-    allocate(hdr%int4tags(size(i_int4)),hdr%int4vals(size(i_int4)),stat=ierrs(4))
+ if (size_type(i_int4) > 0) then
+    allocate(hdr%int4tags(size_type(i_int4)),hdr%int4vals(size_type(i_int4)),stat=ierrs(4))
     if (ierrs(4)==0) hdr%int4tags(:) = ''
  endif
- if (size(i_int8) > 0) then
-    allocate(hdr%int8tags(size(i_int8)),hdr%int8vals(size(i_int8)),stat=ierrs(5))
+ if (size_type(i_int8) > 0) then
+    allocate(hdr%int8tags(size_type(i_int8)),hdr%int8vals(size_type(i_int8)),stat=ierrs(5))
     if (ierrs(5)==0) hdr%int8tags(:) = ''
  endif
- if (size(i_real) > 0)  then
-    allocate(hdr%realtags(size(i_real)),hdr%realvals(size(i_real)),stat=ierrs(6))
+ if (size_type(i_real) > 0)  then
+    allocate(hdr%realtags(size_type(i_real)),hdr%realvals(size_type(i_real)),stat=ierrs(6))
     if (ierrs(6)==0) hdr%realtags(:) = ''
  endif
- if (size(i_real4) > 0)  then
-    allocate(hdr%real4tags(size(i_real4)),hdr%real4vals(size(i_real4)),stat=ierrs(7))
+ if (size_type(i_real4) > 0)  then
+    allocate(hdr%real4tags(size_type(i_real4)),hdr%real4vals(size_type(i_real4)),stat=ierrs(7))
     if (ierrs(7)==0) hdr%real4tags(:) = ''
  endif
- if (size(i_real8) > 0)  then
-    allocate(hdr%real8tags(size(i_real8)),hdr%real8vals(size(i_real8)),stat=ierrs(8))
+ if (size_type(i_real8) > 0)  then
+    allocate(hdr%real8tags(size_type(i_real8)),hdr%real8vals(size_type(i_real8)),stat=ierrs(8))
     if (ierrs(8)==0) hdr%real8tags(:) = ''
  endif
 

src/setup/setup_cluster.f90

@@ -112,9 +112,6 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
     dist_fac    = 1.0     ! distance code unit: dist_fac * pc
  endif
 
- !--Set units
- call set_units(dist=dist_fac*pc,mass=mass_fac*solarm,G=1.)
-
  if (maxvxyzu >= 4) ieos_in = 2 ! Adiabatic equation of state
 
  !--Read values from .setup
@@ -131,6 +128,9 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
     call write_setupfile(fileset)
  endif
 
+ !--Set units
+ call set_units(dist=dist_fac*pc,mass=mass_fac*solarm,G=1.)
+
  !--Define remaining variables using the inputs
  polyk         = kboltz*Temperature/(mu*mass_proton_cgs)*(utime/udist)**2
  rmax          = Rcloud_pc*(pc/udist)

src/tests/test_ptmass.f90

@@ -766,7 +766,7 @@ subroutine test_createsink(ntests,npass)
  use part,       only:init_part,npart,npartoftype,igas,xyzh,massoftype,hfact,rhoh,&
                       iphase,isetphase,fext,divcurlv,vxyzu,fxyzu,poten, &
                       nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,ndptmass, &
-                      dptmass
+                      dptmass,fxyz_ptmass_sinksink
  use ptmass,     only:ptmass_accrete,update_ptmass,icreate_sinks,&
                       ptmass_create,finish_ptmass,ipart_rhomax,h_acc,rho_crit,rho_crit_cgs
  use energies,   only:compute_energies,angtot,etot,totmom
@@ -886,7 +886,7 @@ subroutine test_createsink(ntests,npass)
        call reduceloc_mpi('max',ipart_rhomax_global,id_rhomax)
     endif
     call ptmass_create(nptmass,npart,itestp,xyzh,vxyzu,fxyzu,fext,divcurlv,poten,&
-                       massoftype,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dptmass,0.)
+                       massoftype,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,fxyz_ptmass_sinksink,dptmass,0.)
     !
     ! check that creation succeeded
     !