Mass transfer setup

build/Makefile_setups

@@ -845,6 +845,15 @@ ifeq ($(SETUP), windtunnel)
     ANALYSIS=analysis_common_envelope.f90
 endif
 
+ifeq ($(SETUP), masstransfer)
+#   Wind tunnel setup
+    SETUPFILE= setup_masstransfer.f90
+    GRAVITY=yes
+    KNOWN_SETUP=yes
+    IND_TIMESTEPS=yes
+    ANALYSIS=analysis_common_envelope.f90
+endif
+
 ifeq ($(SETUP), jet)
 #   Jet simulation from Price, Tricco & Bate (2012)
     SETUPFILE= velfield_fromcubes.f90 setup_sphereinbox.f90

src/main/eos.f90

@@ -1444,13 +1444,17 @@ subroutine write_options_eos(iunit)
  endif
 
  select case(ieos)
+ case(2)
+    call write_inopt(gamma,'gamma','Adiabatic index',iunit)
  case(8)
     call write_options_eos_barotropic(iunit)
  case(9)
     call write_options_eos_piecewise(iunit)
  case(10)
     call write_inopt(X_in,'X','hydrogen mass fraction',iunit)
     call write_inopt(Z_in,'Z','metallicity',iunit)
+ case(12)
+    call write_inopt(gamma,'gamma','Adiabatic index',iunit) 
  case(15) ! helmholtz eos
     call eos_helmholtz_write_inopt(iunit)
  case(20)

src/setup/set_binary.f90

@@ -25,7 +25,7 @@ module setbinary
 ! :Dependencies: binaryutils
 !
  implicit none
- public :: set_binary,Rochelobe_estimate,L1_point,get_a_from_period
+ public :: set_binary,Rochelobe_estimate,L1_point,get_a_from_period,get_period_from_a
  public :: get_mean_angmom_vector,get_eccentricity_vector
 
  private
@@ -406,6 +406,17 @@ function get_a_from_period(m1,m2,period) result(a)
 
 end function get_a_from_period
 
+!-------------------------------------------------------------
+! Function to determine the period given the semi-major axis
+!-------------------------------------------------------------
+function get_period_from_a(m1,m2,a) result(period)
+ real, intent(in) :: m1,m2,a
+ real :: period
+
+period= sqrt(((2.*pi)**2*a**3)/(m1 + m2))
+
+end function get_period_from_a
+
 !----------------------------------------------------
 ! Eccentricity vector, for second body
 ! https://en.wikipedia.org/wiki/Eccentricity_vector

src/setup/setup_masstransfer.f90

@@ -0,0 +1,215 @@
+!--------------------------------------------------------------------------!
+! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. !
+! Copyright (c) 2007-2024 The Authors (see AUTHORS)                        !
+! See LICENCE file for usage and distribution conditions                   !
+! http://phantomsph.github.io/                                             !
+!--------------------------------------------------------------------------!
+module setup
+!
+! Setup of two stars or sink particles in a binary
+!
+! :References: None
+!
+! :Owner: Daniel Price
+!
+! :Runtime parameters:
+!   - a            : *semi-major axis*
+!   - mdon         : *donor/primary star mass*
+!   - macc         : *accretor/companion star mass* 
+!   - corotate     : *set stars in corotation*
+!   - eccentricity : *eccentricity*
+!   - f            : *initial true anomaly (180=apoastron)*
+!   - inc          : *inclination (deg)*
+!   - relax        : *relax stars into equilibrium*
+!   - w            : *argument of periapsis (deg)*
+!
+! :Dependencies: centreofmass, dim, eos, externalforces, infile_utils, io,
+!   mpidomain, options, part, physcon, relaxstar, setbinary, setstar,
+!   setunits, setup_params, units
+!
+
+ implicit none
+ public :: setpart
+ real    :: a,mdon,macc,hacc
+
+ private
+
+contains
+
+!----------------------------------------------------------------
+!+
+!  setup for binary star simulations (with or without gas)
+!+
+!----------------------------------------------------------------
+subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,&
+                   polyk,gamma,hfact,time,fileprefix)
+ use part,           only:nptmass,xyzmh_ptmass,vxyz_ptmass,ihacc
+ use setbinary,      only:set_binary,get_period_from_a
+use centreofmass,    only:reset_centreofmass
+ use options,        only:iexternalforce
+ use externalforces, only:iext_corotate,omega_corotate
+ use extern_corotate, only:icompanion_grav,companion_xpos,companion_mass,hsoft
+ use io,             only:master,fatal
+ use eos,            only:ieos, gmw
+ use setunits,       only:mass_unit,dist_unit
+ use timestep,       only:tmax,dtmax
+ integer,           intent(in)    :: id
+ integer,           intent(inout) :: npart
+ integer,           intent(out)   :: npartoftype(:)
+ real,              intent(out)   :: xyzh(:,:)
+ real,              intent(out)   :: massoftype(:)
+ real,              intent(out)   :: polyk,gamma,hfact
+ real,              intent(inout) :: time
+ character(len=20), intent(in)    :: fileprefix
+ real,              intent(out)   :: vxyzu(:,:)
+ character(len=120) :: filename
+ integer :: ierr
+ logical :: iexist
+ real    :: period,ecc,hdon,mass_ratio
+!
+!--general parameters
+!
+ dist_unit = 'solarr'
+ mass_unit = 'solarm'
+ time = 0.
+ polyk = 0.
+ gamma = 5./3.
+!
+!--space available for injected gas particles
+!  in case only sink particles are used
+!
+ npart = 0
+ npartoftype(:) = 0
+ massoftype = 0.
+
+ iexternalforce = iext_corotate
+ icompanion_grav = 1
+ xyzh(:,:)  = 0.
+ vxyzu(:,:) = 0.
+ nptmass = 0
+ a    = 266.34
+ mdon = 6.97
+ macc = 1.41
+ hacc = 1.
+ ieos = 2
+ gmw  = 0.6
+ ecc  = 0.
+ hdon = 1.
+
+ if (id==master) print "(/,65('-'),1(/,a),/,65('-'),/)",&
+   ' Welcome to the Ultimate Binary Setup'
+
+ filename = trim(fileprefix)//'.setup'
+ inquire(file=filename,exist=iexist)
+ if (iexist) call read_setupfile(filename,ieos,polyk,ierr)
+ if (.not. iexist .or. ierr /= 0) then
+    if (id==master) then
+       call write_setupfile(filename)
+       print*,' Edit '//trim(filename)//' and rerun phantomsetup'
+    endif
+    stop
+ endif
+ !
+ !
+ !--if a is negative or is given time units, interpret this as a period
+ !
+
+ period = get_period_from_a(mdon,macc,a)
+ tmax = 10.*period
+ dtmax = tmax/200.
+ !
+ !--now setup orbit using fake sink particles
+ !
+ call set_binary(mdon,macc,a,ecc,hdon,hacc,xyzmh_ptmass,vxyz_ptmass,nptmass,ierr,omega_corotate,&
+                  verbose=(id==master))
+
+ call reset_centreofmass(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass)
+
+
+ if (ierr /= 0) call fatal ('setup_binary','error in call to set_binary')
+
+ companion_mass = mdon
+ companion_xpos = xyzmh_ptmass(1,1)
+ mass_ratio = mdon / macc
+ hsoft = 0.1 * 0.49 * mass_ratio**(2./3.) / (0.6*mass_ratio**(2./3.) + &
+               log( 1. + mass_ratio**(1./3.) ) ) * a
+ !
+ !--delete donor sink 
+ !
+ nptmass=1
+ xyzmh_ptmass(:,1) = xyzmh_ptmass(:,2)
+ vxyz_ptmass(1:3,1) = 0.
+
+ !--restore options
+ !
+
+
+end subroutine setpart
+
+!----------------------------------------------------------------
+!+
+!  write options to .setup file
+!+
+!----------------------------------------------------------------
+subroutine write_setupfile(filename)
+ use infile_utils, only:write_inopt
+ use setunits,     only:write_options_units
+ use eos,          only:write_options_eos
+ character(len=*), intent(in) :: filename
+ integer :: iunit
+
+ print "(a)",' writing setup options file '//trim(filename)
+ open(newunit=iunit,file=filename,status='replace',form='formatted')
+ write(iunit,"(a)") '# input file for binary setup routines'
+
+ call write_options_units(iunit)
+ call write_options_eos(iunit)
+
+ write(iunit,"(/,a)") '# orbit settings'
+ call write_inopt(a,'a','semi-major axis',iunit)
+ call write_inopt(mdon,'mdon','mass of the donor star',iunit)
+ call write_inopt(macc,'macc','mass of the companion star',iunit)
+ call write_inopt(hacc,'hacc','accretion radius of the companion star',iunit)
+
+ close(iunit)
+
+end subroutine write_setupfile
+
+!----------------------------------------------------------------
+!+
+!  read options from .setup file
+!+
+!----------------------------------------------------------------
+subroutine read_setupfile(filename,ieos,polyk,ierr)
+ use infile_utils, only:open_db_from_file,inopts,read_inopt,close_db
+ use io,           only:error,fatal
+ use setunits,     only:read_options_and_set_units
+ character(len=*), intent(in) :: filename
+ integer,          intent(inout) :: ieos
+ real,             intent(inout) :: polyk
+ integer,          intent(out) :: ierr
+ integer, parameter :: iunit = 21
+ integer :: nerr
+ type(inopts), allocatable :: db(:)
+
+ nerr = 0
+ ierr = 0 
+
+ call open_db_from_file(db,filename,iunit,ierr)
+ call read_options_and_set_units(db,nerr)
+
+ call read_inopt(ieos,'ieos',db,errcount=nerr) ! equation of state
+ call read_inopt(a,'a',db,errcount=nerr)
+ call read_inopt(mdon,'mdon',db,errcount=nerr)
+ call read_inopt(macc,'macc',db,errcount=nerr)
+ call read_inopt(hacc,'hacc',db,errcount=nerr)
+ call close_db(db)
+
+  if (nerr > 0) then
+    print "(1x,i2,a)",nerr,' error(s) during read of setup file: re-writing...'
+    ierr = nerr
+ endif
+
+end subroutine read_setupfile
+
+end module setup