Boost factor for gravitational-wave inspiral

src/main/extern_gwinspiral.f90

@@ -24,11 +24,12 @@ module extern_gwinspiral
  ! Runtime parameters
  !
  real,    public  :: stopratio = 0.005
+ real,    public  :: gw_boostfac = 1.  ! multiplier to force applied to stars
  !
  ! local variables
  !
  integer, private :: n_threshhold
- integer, public  :: Nstar(2) = 0  ! give default value in case dump header not read
+ integer, public  :: Nstar_gw(2) = 0  ! give default value in case dump header not read
  real,    private :: fstar1_coef,fstar2_coef
  real,    private :: com(3),comstar1(3),comstar2(3),vcomstar1(3),vcomstar2(3),fstar1(3),fstar2(3)
  logical, private :: isseparate = .true.
@@ -60,12 +61,12 @@ subroutine initialise_gwinspiral(npart,nptmass,ierr)
 
  ierr = 0
  nerr = 0
- if (Nstar(1) > 0) then
-    write(*,"(2(a,i8))") ' Initialising inspiral two stars scenario, Nstar_1 = ',Nstar(1),' Nstar_2 = ', Nstar(2)
- elseif (Nstar(1)==0 .and. Nstar(2)==0 .and. nptmass==2) then
+ if (Nstar_gw(1) > 0) then
+    write(*,"(2(a,i8))") ' Initialising inspiral two stars scenario, Nstar_1 = ',Nstar_gw(1),' Nstar_2 = ', Nstar_gw(2)
+ elseif (Nstar_gw(1)==0 .and. Nstar_gw(2)==0 .and. nptmass==2) then
     write(*,"(2(a,i8))") ' Initialising inspiral on two sink particles'
- elseif (Nstar(2)==0 .and. nptmass==1) then
-    write(*,"(2(a,i8))") ' Initialising inspiral star-particle scenario, Nstar_1 = ',Nstar(1)
+ elseif (Nstar_gw(2)==0 .and. nptmass==1) then
+    write(*,"(2(a,i8))") ' Initialising inspiral star-particle scenario, Nstar_1 = ',Nstar_gw(1)
  else
     ierr = 1
     isseparate = .false.
@@ -100,21 +101,21 @@ subroutine gw_still_inspiralling(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptma
  fstar2_coef = 0.
  stopped_now = .false.
  if ( isseparate ) then
-    if (Nstar(1) == 0 .and. nptmass==2) then
+    if (Nstar_gw(1) == 0 .and. nptmass==2) then
        comstar1  = xyzmh_ptmass(1:3,2)
        mstar1    = xyzmh_ptmass(4,2)
        vcomstar1 = vxyz_ptmass(1:3,2)
     else
-       call get_centreofmass(comstar1,vcomstar1,Nstar(1),xyzh(:,1:Nstar(1)),&
-                          vxyzu(:,1:Nstar(1)),mass=mstar1)
+       call get_centreofmass(comstar1,vcomstar1,Nstar_gw(1),xyzh(:,1:Nstar_gw(1)),&
+                          vxyzu(:,1:Nstar_gw(1)),mass=mstar1)
     endif
-    if (Nstar(2) == 0 .and. nptmass>=1) then
+    if (Nstar_gw(2) == 0 .and. nptmass>=1) then
        comstar2  = xyzmh_ptmass(1:3,1)
        mstar2    = xyzmh_ptmass(4,1)
        vcomstar2 = vxyz_ptmass(1:3,1)
     else
-       call get_centreofmass(comstar2,vcomstar2,Nstar(2),xyzh(:,Nstar(1)+1:npart),&
-                          vxyzu(:,Nstar(1)+1:npart),mass=mstar2)
+       call get_centreofmass(comstar2,vcomstar2,Nstar_gw(2),xyzh(:,Nstar_gw(1)+1:npart),&
+                          vxyzu(:,Nstar_gw(1)+1:npart),mass=mstar2)
     endif
     com  = (comstar1*mstar1 + comstar2*mstar2)/(mstar1+mstar2)
     dirstar1 = comstar1 - com ! The directional vector to star 1
@@ -124,11 +125,11 @@ subroutine gw_still_inspiralling(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptma
     !
     k1 = 0
     !$omp parallel default(none) &
-    !$omp shared(nstar,xyzh,com,dirstar1) &
+    !$omp shared(Nstar_gw,xyzh,com,dirstar1) &
     !$omp private(i,dx,dy,dz,dir) &
     !$omp reduction(+:k1)
     !$omp do
-    do i=1,nstar(1)
+    do i=1,Nstar_gw(1)
        dx  = xyzh(1,i) - com(1)
        dy  = xyzh(2,i) - com(2)
        dz  = xyzh(3,i) - com(3)
@@ -143,11 +144,11 @@ subroutine gw_still_inspiralling(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptma
     !
     k2 = 0
     !$omp parallel default(none) &
-    !$omp shared(nstar,npart,xyzh,com,dirstar1) &
+    !$omp shared(Nstar_gw,npart,xyzh,com,dirstar1) &
     !$omp private(i,dx,dy,dz,dir) &
     !$omp reduction(+:k2)
     !$omp do
-    do i=nstar(1)+1,nstar(1)+nstar(2)
+    do i=Nstar_gw(1)+1,Nstar_gw(1)+Nstar_gw(2)
        dx  = xyzh(1,i) - com(1)
        dy  = xyzh(2,i) - com(2)
        dz  = xyzh(3,i) - com(3)
@@ -203,8 +204,8 @@ subroutine get_gw_force()
     !
     ! Compute the drag force vectors for each star
     !
-    fstar1 = fstar1_coef * vcomstar1 / (vstar1sq*separation**5)
-    fstar2 = fstar2_coef * vcomstar2 / (vstar2sq*separation**5)
+    fstar1 = gw_boostfac * fstar1_coef * vcomstar1 / (vstar1sq*separation**5)
+    fstar2 = gw_boostfac * fstar2_coef * vcomstar2 / (vstar2sq*separation**5)
  endif
 
 end subroutine get_gw_force
@@ -220,20 +221,20 @@ subroutine get_gw_force_i(i,fextxi,fextyi,fextzi,phi)
  real,    intent(inout) :: fextxi,fextyi,fextzi,phi
 
  if (i > 0 .and. isseparate ) then
-    if (i <= nstar(1)) then
+    if (i <= Nstar_gw(1)) then
        fextxi = fstar1(1)
        fextyi = fstar1(2)
        fextzi = fstar1(3)
-    elseif (nstar(2) > 0) then
+    elseif (Nstar_gw(2) > 0) then
        fextxi = fstar2(1)
        fextyi = fstar2(2)
        fextzi = fstar2(3)
     endif
- elseif (i == -1 .and. nstar(2)==0 .and. isseparate) then
+ elseif (i == -1 .and. Nstar_gw(2)==0 .and. isseparate) then
     fextxi = fstar2(1)  ! acceleration applied to sink particle 1 (star 2)
     fextyi = fstar2(2)
     fextzi = fstar2(3)
- elseif (i == -2 .and. nstar(1)==0 .and. isseparate) then
+ elseif (i == -2 .and. Nstar_gw(1)==0 .and. isseparate) then
     fextxi = fstar1(1) ! acceleration applied to sink particle 2 (star 1)
     fextyi = fstar1(2)
     fextzi = fstar1(3)
@@ -255,6 +256,7 @@ subroutine write_options_gwinspiral(iunit)
  integer, intent(in) :: iunit
 
  call write_inopt(stopratio,'stop_ratio','ratio of particles crossing CoM to indicate a merger',iunit)
+ call write_inopt(gw_boostfac,'gw_boostfac','multiplicative factor to boost GW force',iunit)
 
 end subroutine write_options_gwinspiral
 
@@ -269,6 +271,7 @@ subroutine read_options_gwinspiral(db,nerr)
  integer,      intent(inout) :: nerr
 
  call read_inopt(stopratio,'stop_ratio',db,errcount=nerr,min=0.,max=1.)
+ call read_inopt(gw_boostfac,'gw_boostfac',db,errcount=nerr,min=0.)
 
 end subroutine read_options_gwinspiral
 
@@ -283,8 +286,8 @@ subroutine write_headeropts_gwinspiral(hdr,ierr)
  integer,      intent(out)   :: ierr
 
  ierr = 0
- call add_to_header(Nstar(1),'Nstar_1',hdr,ierr)
- call add_to_header(Nstar(2),'Nstar_2',hdr,ierr)
+ call add_to_header(Nstar_gw(1),'Nstar_1',hdr,ierr)
+ call add_to_header(Nstar_gw(2),'Nstar_2',hdr,ierr)
 
 end subroutine write_headeropts_gwinspiral
 
@@ -301,8 +304,8 @@ subroutine read_headeropts_gwinspiral(hdr,nptmass,ierr)
  integer :: ierr1,ierr2
 
  ierr  = 0
- call extract('Nstar_1',Nstar(1),hdr,ierr1)
- call extract('Nstar_2',Nstar(2),hdr,ierr2)
+ call extract('Nstar_1',Nstar_gw(1),hdr,ierr1)
+ call extract('Nstar_2',Nstar_gw(2),hdr,ierr2)
 
  if (ierr1 /= 0 .or. ierr2 /= 0) then
     ! if there are two sink particles and Nstar_1

src/main/externalforces.f90

@@ -818,7 +818,8 @@ subroutine initialise_externalforces(iexternalforce,ierr)
  end select
 
  select case(iexternalforce)
- case(iext_star,iext_binary,iext_corot_binary,iext_prdrag,iext_spiral,iext_lensethirring,iext_einsteinprec,iext_gnewton)
+ case(iext_star,iext_binary,iext_corot_binary,iext_prdrag,iext_spiral,iext_lensethirring,&
+      iext_einsteinprec,iext_gnewton,iext_gwinspiral)
     !
     !--check that G=1 in code units
     !

src/setup/setup_binary.f90

@@ -30,7 +30,7 @@ module setup
  implicit none
  public :: setpart
 
- logical :: relax,write_rho_to_file,corotate
+ logical :: relax,write_rho_to_file,corotate,gwinspiral
  integer, parameter :: max_stars = max_hier_levels
  type(star_t)  :: star(max_stars)
  type(orbit_t) :: orbit
@@ -53,7 +53,8 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,&
  use setorbit,       only:set_defaults_orbit,set_orbit
  use sethierarchical,only:set_hierarchical,set_hierarchical_default_options,hs,generate_hierarchy_string
  use options,        only:iexternalforce,alphau
- use externalforces, only:iext_corotate,iext_geopot,iext_star,omega_corotate,mass1,accradius1
+ use externalforces, only:iext_corotate,iext_geopot,iext_star,iext_gwinspiral,omega_corotate,mass1,accradius1
+ use extern_gwinspiral, only:Nstar_gw
  use io,             only:master,fatal
  use setstar,        only:set_defaults_stars,set_stars,shift_stars
  use eos,            only:X_in,Z_in,use_var_comp
@@ -109,6 +110,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,&
  call set_hierarchical_default_options()
  relax = .true.
  corotate = .false.
+ gwinspiral = .false.
  use_var_comp = .false.
  write_rho_to_file = .true.
 
@@ -181,6 +183,10 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,&
     accradius1 = xyzmh_ptmass(ihacc,nptmass+1)
     xyzmh_ptmass(:,nptmass) = xyzmh_ptmass(:,nptmass+1)
     vxyz_ptmass(:,nptmass) = vxyz_ptmass(:,nptmass+1)
+ elseif (nstar == 2 .and. (gwinspiral .or. iexternalforce==iext_gwinspiral)) then
+    Nstar_gw(1) = star(1)%np
+    Nstar_gw(2) = npart - star(1)%np
+    iexternalforce = iext_gwinspiral
  elseif (nstar >= 1 .and. set_oblateness .and. star(1)%iprofile == 0) then
     ! set J2 for sink particle 1 to be equal to oblateness of Saturn
     xyzmh_ptmass(iJ2,1) = 0.01629
@@ -222,6 +228,7 @@ subroutine write_setupfile(filename)
     call write_hierarchical_setupfile(iunit,nstar)
  elseif (nstar == 2) then
     call write_inopt(corotate,'corotate','set stars in corotation',iunit)
+    call write_inopt(gwinspiral,'gwinspiral','set gravitational wave inspiral',iunit)
     call write_options_orbit(orbit,iunit)
 
     write(iunit,"(/,a)") '# timestepping'
@@ -261,6 +268,7 @@ subroutine read_setupfile(filename,ierr)
     call read_hierarchical_setupfile(db,nerr,nstar)
  elseif (nstar == 2) then
     call read_inopt(corotate,'corotate',db,errcount=nerr)
+    call read_inopt(gwinspiral,'gwinspiral',db,errcount=nerr)
     m1 = in_code_units(star(1)%m,ierr,unit_type='mass')
     m2 = in_code_units(star(2)%m,ierr1,unit_type='mass')
     if (ierr /= 0 .or. ierr1 /= 0) then

src/utils/analysis_NSmerger.f90

@@ -25,7 +25,7 @@ module analysis
  use physcon,         only: pi
  use centreofmass,    only: get_centreofmass
  use readwrite_dumps, only:opened_full_dump
- use extern_gwinspiral, only:Nstar
+ use extern_gwinspiral, only:Nstar_gw
  implicit none
  character(len=20), parameter, public :: analysistype = 'NSmerger'
  !
@@ -75,8 +75,8 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit)
     !
     !--Get the index range for each star
     ! (Note from DJP: This should now be automatically read from the dump header)
-    if (Nstar(1) <= 0) call fatal('analysis_NSmerger','Require Nstar(1) > 0 in header of dump file')
-    if (Nstar(2) <= 0) call fatal('analysis_NSmerger','Require Nstar(2) > 0 in header of dump file')
+    if (Nstar_gw(1) <= 0) call fatal('analysis_NSmerger','Require Nstar_gw(1) > 0 in header of dump file')
+    if (Nstar_gw(2) <= 0) call fatal('analysis_NSmerger','Require Nstar_gw(2) > 0 in header of dump file')
     !
     !--Prompt for the density_cut off
     if (choice > 1) then
@@ -148,8 +148,8 @@ subroutine trace_com(dumpfile,xyzh,vxyzu,time,npart,iunit)
  endif
  !
  !--Get centre of masses of the stars
- call get_centreofmass(com1,vcom1,nstar(1),xyzh(:,1:nstar(1)),vxyzu(:,1:nstar(1)))
- call get_centreofmass(com2,vcom2,nstar(2),xyzh(:,nstar(1)+1:npart),vxyzu(:,nstar(1)+1:npart))
+ call get_centreofmass(com1,vcom1,Nstar_gw(1),xyzh(:,1:Nstar_gw(1)),vxyzu(:,1:Nstar_gw(1)))
+ call get_centreofmass(com2,vcom2,Nstar_gw(2),xyzh(:,Nstar_gw(1)+1:npart),vxyzu(:,Nstar_gw(1)+1:npart))
  !
  rad = sqrt( (com1(1)-com2(1))**2 + (com1(2)-com2(2))**2 + (com1(3)-com2(3))**2 )
  !

src/utils/moddump_binary.f90

@@ -35,9 +35,10 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
  use centreofmass,      only:reset_centreofmass,get_centreofmass
  use prompting,         only:prompt
  use options,           only:iexternalforce
- use externalforces,    only:omega_corotate,iext_corotate
+ use externalforces,    only:omega_corotate,iext_corotate,iext_gwinspiral
  use extern_corotate,   only:icompanion_grav,companion_xpos,companion_mass,primarycore_xpos,&
                              primarycore_mass,primarycore_hsoft,hsoft
+ use extern_gwinspiral, only:Nstar_gw
  use infile_utils,      only:open_db_from_file,inopts,read_inopt,close_db
  use table_utils,       only:yinterp
  use readwrite_mesa,    only:read_mesa
@@ -62,7 +63,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
  real                      :: time2,hfact2
  real                      :: xyzmh1_stash(nsinkproperties),xyzmh2_stash(nsinkproperties),vxyz1_stash(3),vxyz2_stash(3)
  real, allocatable         :: r(:),den(:),pres(:),temp(:),enitab(:),Xfrac(:),Yfrac(:),mu(:),m(:)
- logical                   :: use_corotating_frame,iprimary_grav_ans
+ logical                   :: use_corotating_frame,iprimary_grav_ans,gwinspiral
  character(len=20)         :: filename = 'binary.in'
  character(len=100)        :: densityfile,dumpname
  type(inopts), allocatable :: db(:)
@@ -178,6 +179,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
        hacc = 0.
        hsoft = 0.
        use_corotating_frame = .false.
+       gwinspiral = .false.
 
        call reset_centreofmass(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass)
        call delete_dead_or_accreted_particles(npart,npartoftype)  !removes the dead or accreted particles for a correct total mass computation
@@ -206,6 +208,9 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
           call prompt('Enter softening length for companion', hsoft, 0.)
        endif
        call prompt('Do you want to transform to a corotating frame and simulate corotating binary?', use_corotating_frame)
+       call prompt('Do you want to add gravitational radiation reaction?', gwinspiral)
+
+       if (gwinspiral) iexternalforce = iext_gwinspiral
 
        ! set the binary
        if (use_corotating_frame) then
@@ -302,13 +307,20 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
              vxyz_ptmass(1:3,nptmass1+nptmass2) = vxyz2_stash(1:3)
           endif
 
+          if (gwinspiral) then
+             Nstar_gw(1) = nstar1
+             Nstar_gw(2) = nstar2
+          endif
+
        else
           nptmass = nptmass1 + 1
           xyzmh_ptmass(1:3,nptmass) = xyzmh2_stash(1:3)
           vxyz_ptmass(1:3,nptmass) = vxyz2_stash(1:3)
           xyzmh_ptmass(4,nptmass) = m2
           xyzmh_ptmass(ihacc,nptmass) = hacc
           xyzmh_ptmass(ihsoft,nptmass) = hsoft
+
+          if (gwinspiral) Nstar_gw(1) = nstar1
        endif
 
        if (nptmass1 == 1) then

src/utils/moddump_sinkbinary.f90

@@ -32,7 +32,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
  use timestep, only:dtmax,tmax
  use options,   only:iexternalforce
  use externalforces, only:iext_gwinspiral
- use extern_gwinspiral, only:Nstar
+ use extern_gwinspiral, only:Nstar_gw
  integer, intent(inout) :: npart
  integer, intent(inout) :: npartoftype(:)
  real,    intent(inout) :: massoftype(:)
@@ -77,7 +77,7 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
  call reset_centreofmass(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass)
 
  if (iexternalforce==iext_gwinspiral) then
-    Nstar(1) = npart
+    Nstar_gw(1) = npart
  endif
 
  period = 2.*pi*sqrt(a**3/(m1 + m2))