Merge pull request #583 from Yrisch/Slow-down

Slow down method and few fixes on the new star formation prescription

.mailmap

@@ -123,5 +123,5 @@ Nicolás Cuello <[email protected]> Nicolas Cuello <[email protected]
 Rebecca Martin <[email protected]> rebeccagmartin <[email protected]>
 Stephen Nielson <[email protected]> s-neilson <[email protected]>
 Stephen Nielson <[email protected]> Stephen Neilson <[email protected]>
-Yann Bernard <[email protected]> Yrisch <[email protected]>
+Yann Bernard <[email protected]> Yrisch <[email protected]>
 David Bamba <[email protected]> DavidBamba <[email protected]>

AUTHORS

@@ -12,8 +12,8 @@ James Wurster <[email protected]>
 David Liptai <[email protected]>
 Lionel Siess <[email protected]>
 Fangyi (Fitz) Hu <[email protected]>
+Yann Bernard <[email protected]>
 Daniel Mentiplay <[email protected]>
-Yann Bernard <[email protected]>
 Megha Sharma <[email protected]>
 Arnaud Vericel <[email protected]>
 Mark Hutchison <[email protected]>
@@ -42,35 +42,36 @@ Sahl Rowther <[email protected]>
 Simon Glover <[email protected]>
 Thomas Reichardt <[email protected]>
 Jean-François Gonzalez <[email protected]>
-Madeline Overton <[email protected]>
 Christopher Russell <[email protected]>
-Jolien Malfait <[email protected]>
-Phantom benchmark bot <[email protected]>
+Madeline Overton <[email protected]>
 Alessia Franchini <[email protected]>
 Alex Pettitt <[email protected]>
+Jolien Malfait <[email protected]>
+Phantom benchmark bot <[email protected]>
 Kieran Hirsh <[email protected]>
 Nicole Rodrigues <[email protected]>
+Ana Lourdes Juarez <[email protected]>
+David Trevascus <[email protected]>
 Farzana Meru <[email protected]>
 Nicolás Cuello <[email protected]>
-David Trevascus <[email protected]>
-Miguel Gonzalez-Bolivar <[email protected]>
 Chris Nixon <[email protected]>
-Shunquan Huang <[email protected]>
-Joe Fisher <[email protected]>
+Miguel Gonzalez-Bolivar <[email protected]>
 Benoit Commercon <[email protected]>
-Zachary Pellow <[email protected]>
 Giulia Ballabio <[email protected]>
+Joe Fisher <[email protected]>
 Maxime Lombart <[email protected]>
 Orsola De Marco <[email protected]>
-Steven Rieder <[email protected]>
-Stéven Toupin <[email protected]>
-Taj Jankovič <[email protected]>
-Jeremy Smallwood <[email protected]>
+Shunquan Huang <[email protected]>
+Zachary Pellow <[email protected]>
 Ariel Chitan <[email protected]>
-Rebecca Martin <[email protected]>
-Jorge Cuadra <[email protected]>
-David Bamba <[email protected]>
-Shunquan Huang <[email protected]>
-Cox, Samuel <[email protected]>
 Chunliang Mu <[email protected]>
+Cox, Samuel <[email protected]>
+David Bamba <[email protected]>
 Hugh Griffiths <[email protected]>
+Jeremy Smallwood <[email protected]>
+Jorge Cuadra <[email protected]>
+Rebecca Martin <[email protected]>
+Shunquan Huang <[email protected]>
+Steven Rieder <[email protected]>
+Stéven Toupin <[email protected]>
+Taj Jankovič <[email protected]>

build/Makefile

@@ -525,7 +525,7 @@ SOURCES= physcon.f90 ${CONFIG} ${SRCKERNEL} io.F90 units.f90 \
          utils_system.f90 utils_mathfunc.f90 part.F90 ${DOMAIN} boundary.f90 boundary_dynamic.f90 utils_timing.f90 mpi_balance.F90 \
          setup_params.f90 timestep.f90 utils_dumpfiles.f90 utils_indtimesteps.F90 \
          utils_sort.f90 utils_supertimestep.F90 utils_tables.f90 utils_gravwave.f90 \
-         utils_sphNG.f90 utils_vectors.f90 utils_datafiles.f90 datafiles.f90 \
+         utils_sphNG.f90 utils_vectors.f90 utils_subgroup.f90 utils_kepler.f90 utils_datafiles.f90 datafiles.f90 \
          gitinfo.f90 ${SRCFASTMATH} random.f90 ${SRCEOS} cullendehnen.f90 ${SRCNIMHD} ${SRCGR} \
          checkoptions.F90 viscosity.f90 damping.f90 options.f90 cons2primsolver.f90 radiation_utils.f90 cons2prim.f90 \
          centreofmass.f90 ${SRCPOT} checkconserved.f90 \
@@ -535,7 +535,7 @@ SOURCES= physcon.f90 ${CONFIG} ${SRCKERNEL} io.F90 units.f90 \
          utils_deriv.f90 utils_implicit.f90 radiation_implicit.f90 ${SRCTURB} \
          ${SRCINJECT_DEPS} wind_equations.f90 wind.F90  \
          ${SRCKROME} memory.f90 ${SRCREADWRITE_DUMPS} ${SRCINJECT} \
-         H2regions.f90 utils_subgroup.f90 utils_kepler.f90 subgroup.f90 \
+         H2regions.f90 subgroup.f90 \
          quitdump.f90 ptmass.F90\
          readwrite_infile.F90 dens.F90 force.F90 deriv.F90 energies.F90 sort_particles.f90 \
          utils_shuffleparticles.F90 evwrite.f90 substepping.F90 step_leapfrog.F90 writeheader.F90 ${SRCAN} step_supertimestep.F90 \

src/main/H2regions.f90

@@ -6,7 +6,10 @@
 !--------------------------------------------------------------------------!
 module HIIRegion
 !
-! Feedback from HII regions
+! HIIRegion
+! contains routines to model HII region expansion due to ionization and radiation pressure..
+! routine originally made by Hopkins et al. (2012),reused by Fujii et al. (2021)
+! and adapted in Phantom by Yann Bernard
 !
 ! :References: Fujii et al. (2021), Hopkins et al. (2012)
 !

src/main/cons2prim.f90

@@ -20,7 +20,7 @@ module cons2prim
 !   Liptai & Price (2019), MNRAS 485, 819-842
 !   Ballabio et al. (2018), MNRAS 477, 2766-2771
 !
-! :Owner: Elisabeth Borchert
+! :Owner: Daniel Price
 !
 ! :Runtime parameters: None
 !

src/main/energies.F90

@@ -72,7 +72,7 @@ subroutine compute_energies(t)
                           isdead_or_accreted,epot_sinksink,imacc,ispinx,ispiny,&
                           ispinz,mhd,gravity,poten,dustfrac,eos_vars,itemp,igasP,ics,&
                           nden_nimhd,eta_nimhd,iion,ndustsmall,graindens,grainsize,&
-                          iamdust,ndusttypes,rad,iradxi,gtgrad,group_info,n_group
+                          iamdust,ndusttypes,rad,iradxi,gtgrad,group_info,bin_info,n_group
  use part,           only:pxyzu,fxyzu,fext
  use gravwaveutils,  only:calculate_strain,calc_gravitwaves
  use centreofmass,   only:get_centreofmass_accel
@@ -83,7 +83,7 @@ subroutine compute_energies(t)
  use options,        only:iexternalforce,calc_erot,alpha,ieos,use_dustfrac
  use mpiutils,       only:reduceall_mpi
  use ptmass,         only:get_accel_sink_gas,use_regnbody
- use subgroup,     only:get_pot_subsys
+ use subgroup,       only:get_pot_subsys
  use viscosity,      only:irealvisc,shearfunc
  use nicil,          only:nicil_update_nimhd,nicil_get_halldrift,nicil_get_ambidrift, &
                      use_ohm,use_hall,use_ambi,n_data_out,n_warn,eta_constant
@@ -644,7 +644,7 @@ subroutine compute_energies(t)
  erad = reduceall_mpi('+',erad)
  if (nptmass > 1) then
     if (use_regnbody) then
-       call get_pot_subsys(n_group,group_info,xyzmh_ptmass,fxyz_ptmass,gtgrad,epot_sinksink)
+       call get_pot_subsys(n_group,group_info,bin_info,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,gtgrad,epot_sinksink)
     endif
     epot = epot + epot_sinksink
  endif

src/main/evolve.F90

@@ -91,7 +91,7 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag)
                             xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dptmass,gravity,iboundary, &
                             fxyz_ptmass_sinksink,ntot,poten,ndustsmall,accrete_particles_outside_sphere,&
                             linklist_ptmass,isionised,dsdt_ptmass,isdead_or_accreted
- use part,             only:n_group,n_ingroup,n_sing,group_info,nmatrix
+ use part,             only:n_group,n_ingroup,n_sing,group_info,bin_info,nmatrix
  use quitdump,         only:quit
  use ptmass,           only:icreate_sinks,ptmass_create,ipart_rhomax,pt_write_sinkev,calculate_mdot, &
                             set_integration_precision,ptmass_create_stars,use_regnbody,ptmass_create_seeds,&
@@ -143,7 +143,7 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag)
  logical         :: use_global_dt
  integer         :: j,nskip,nskipped,nevwrite_threshold,nskipped_sink,nsinkwrite_threshold
  character(len=120) :: dumpfile_orig
- integer         :: dummy,istepHII
+ integer         :: dummy,istepHII,nptmass_old
 
  dummy = 0
 
@@ -279,7 +279,7 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag)
     !  across all nodes
     nskip = int(ntot)
 #endif
-
+    nptmass_old = nptmass
     if (gravity .and. icreate_sinks > 0 .and. ipart_rhomax /= 0) then
        !
        ! creation of new sink particles
@@ -310,20 +310,21 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag)
           istepHII = 2**nbinmax/HIIuprate
           if (istepHII==0) istepHII = 1
        endif
-       if (mod(istepfrac,istepHII)==0 .or. istepfrac==1 .or. (icreate_sinks == 2 .and. ipart_createstars /= 0)) then
+       if (mod(istepfrac,istepHII) == 0 .or. istepfrac == 1 .or. (icreate_sinks == 2 .and. ipart_createstars /= 0)) then
           call HII_feedback(nptmass,npart,xyzh,xyzmh_ptmass,vxyzu,isionised)
        endif
     endif
 
     ! Need to recompute the force when sink or stars are created
-    if (ipart_rhomax /= 0 .or. ipart_createseeds /= 0 .or. ipart_createstars /= 0) then
+    if (nptmass > nptmass_old .or. ipart_createseeds /= 0 .or. ipart_createstars /= 0) then
        if (use_regnbody) then
-          call group_identify(nptmass,n_group,n_ingroup,n_sing,xyzmh_ptmass,vxyz_ptmass,group_info,nmatrix)
+          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,group_info=group_info)
+                 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)
+                fxyz_ptmass,dsdt_ptmass,0.,0.,dummy,.false.,linklist_ptmass,bin_info)
        endif
        if (ipart_createseeds /= 0) ipart_createseeds = 0 ! reset pointer to zero
        if (ipart_createstars /= 0) ipart_createstars = 0 ! reset pointer to zero

src/main/initial.F90

@@ -140,7 +140,7 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
                             epot_sinksink,get_ntypes,isdead_or_accreted,dustfrac,ddustevol,&
                             nden_nimhd,dustevol,rhoh,gradh, &
                             Bevol,Bxyz,dustprop,filfac,ddustprop,ndustsmall,iboundary,eos_vars,dvdx, &
-                            n_group,n_ingroup,n_sing,nmatrix,group_info,isionised
+                            n_group,n_ingroup,n_sing,nmatrix,group_info,bin_info,isionised
  use part,             only:pxyzu,dens,metrics,rad,radprop,drad,ithick
  use densityforce,     only:densityiterate
  use linklist,         only:set_linklist
@@ -221,7 +221,7 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
  use checkconserved,   only:get_conserv,etot_in,angtot_in,totmom_in,mdust_in
  use fileutils,        only:make_tags_unique
  use damping,          only:idamp
- use subgroup,         only:group_identify,init_subgroup
+ use subgroup,         only:group_identify,init_subgroup,update_kappa
  use HIIRegion,        only:iH2R,initialize_H2R,update_ionrates
  character(len=*), intent(in)  :: infile
  character(len=*), intent(out) :: logfile,evfile,dumpfile
@@ -517,9 +517,11 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
     ! compute initial sink-sink forces and get timestep
     if (use_regnbody) then
        call init_subgroup
-       call group_identify(nptmass,n_group,n_ingroup,n_sing,xyzmh_ptmass,vxyz_ptmass,group_info,nmatrix)
+       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)
+                             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)
@@ -536,8 +538,14 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
           if (ntypes > 1 .and. maxphase==maxp) then
              pmassi = massoftype(iamtype(iphase(i)))
           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)
+          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
           dtsinkgas = min(dtsinkgas,C_force*1./sqrt(fonrmax),C_force*sqrt(dtphi2))
        endif
     enddo
@@ -552,6 +560,7 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
     !  Reduce dt over MPI tasks
     dtsinkgas = reduceall_mpi('min',dtsinkgas)
     dtextforce = reduceall_mpi('min',dtextforce)
+    if (use_regnbody) call update_kappa(xyzmh_ptmass,vxyz_ptmass,bin_info,group_info,n_group)
  endif
  call init_ptmass(nptmass,logfile)
  if (gravity .and. icreate_sinks > 0 .and. id==master) then

src/main/inject_randomwind.f90

@@ -12,7 +12,7 @@ module inject
 ! :References:
 !   Trevascus et al. (2021), MNRAS 505, L21-L25
 !
-! :Owner: Shunquan Huang
+! :Owner: Daniel Price
 !
 ! :Runtime parameters:
 !   - delta_theta : *standard deviation for the gaussion distribution*

src/main/part.F90

@@ -298,11 +298,23 @@ module part
 !
 !-- Regularisation algorithm allocation
 !
- integer,         allocatable :: group_info(:,:)
- integer(kind=1), allocatable :: nmatrix(:,:)
- integer, parameter :: igarg  = 1 ! idx of the particle member of a group
- integer, parameter :: igcum  = 2 ! cumulative sum of the indices to find the starting and ending point of a group
- integer, parameter :: igid   = 3 ! id of the group, correspond to the root of the group in the dfs/union find construction
+ integer(kind=1), allocatable :: nmatrix(:,:)    ! adjacency matrix used to construct each groups
+
+ integer, allocatable :: group_info(:,:) ! array storing group id/idx of each group comp/boundary idx/bin comp id
+ integer, parameter   :: igarg  = 1 ! idx of the particle member of a group
+ integer, parameter   :: igcum  = 2 ! cumulative sum of the indices to find the starting and ending point of a group
+ integer, parameter   :: igid   = 3 ! id of the group, correspond to the root of the group in the dfs/union find construction
+ integer, parameter   :: icomp  = 4 ! id of the binary companion if it exists, otherwise equal to the id
+
+ real, allocatable    :: bin_info(:,:) ! array storing important orbital parameters and quantities of each binary
+ integer, parameter   :: isemi = 1 ! semi major axis
+ integer, parameter   :: iecc  = 2 ! eccentricity
+ integer, parameter   :: iapo  = 3 ! apocenter
+ integer, parameter   :: iorb  = 4 ! orbital period
+ integer, parameter   :: ipert = 5 ! perturbation
+ integer, parameter   :: ikap  = 6 ! kappa slow down
+
+
  ! needed for group identification and sorting
  integer  :: n_group = 0
  integer  :: n_ingroup = 0
@@ -504,7 +516,8 @@ subroutine allocate_part
     call allocate_array('abundance', abundance, nabundances, maxp_h2)
  endif
  call allocate_array('T_gas_cool', T_gas_cool, maxp_krome)
- call allocate_array('group_info', group_info, 3, maxptmass)
+ call allocate_array('group_info', group_info, 4, maxptmass)
+ call allocate_array('bin_info', bin_info, 6, maxptmass)
  call allocate_array("nmatrix", nmatrix, maxptmass, maxptmass)
  call allocate_array("gtgrad", gtgrad, 3, maxptmass)
  call allocate_array('isionised', isionised, maxp)
@@ -590,6 +603,7 @@ subroutine deallocate_part
  if (allocated(istsactive))   deallocate(istsactive)
  if (allocated(ibin_sts))     deallocate(ibin_sts)
  if (allocated(group_info))   deallocate(group_info)
+ if (allocated(bin_info))     deallocate(bin_info)
  if (allocated(nmatrix))      deallocate(nmatrix)
  if (allocated(gtgrad))       deallocate(gtgrad)
  if (allocated(isionised))    deallocate(isionised)