Skip to content

Commit

Permalink
Merge pull request #417 from danieljprice/set_star
Browse files Browse the repository at this point in the history 
restartable relaxation process for stars #414
  • Loading branch information
@danieljprice
danieljprice committed May 12, 2023
2 parents eddddb1 + 8a8024b commit 78e4948
Show file tree
Hide file tree
Showing 40 changed files with 638 additions and 524 deletions.
142 changes: 108 additions & 34 deletions docs/binary.rst
View file
Original file line number Diff line number Diff line change
@@ -1,55 +1,86 @@
Binary stars and common envelope evolution
============================================

Using SETUP=binary
------------------
The one-stop-shop to setup a binary star simulation is to use SETUP=binary::
Setting up and relaxing binary stars (the easy way)
---------------------------------------------------
The one-stop-shop to setup a binary star simulation is as follows::

~/phantom/scripts/writemake.sh binary > Makefile
make setup
./phantomsetup sim
make a new directory and write a local Makefile
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

giving::
::

-----------------------------------------------------------------
$ mkdir mybinarysim
$ cd mybinarysim
$ ~/phantom/scripts/writemake.sh binary > Makefile

compile phantom and phantomsetup
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

::

$ make
$ make setup
$ ls
Makefile phantom* phantomsetup*

run phantomsetup
~~~~~~~~~~~~~~~~

::

./phantomsetup sim

-----------------------------------------------------------------
Welcome to the Ultimate Binary Setup
-----------------------------------------------------------------

writing setup options file binary.setup
Edit sim.setup and rerun phantomsetup


This will create a file called sim.setup which contains setup options. Open this file in
your favourite text editor to proceed...
This will create a file called sim.setup which contains setup options.
Open this file in your favourite text editor to proceed...


Two sink particles in orbit
~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The simplest is to setup two sink particles in a binary (iprofile=0), by amending the iprofile flags::
----------------------------

Set iprofile=0 for both stars in the .setup file
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

::

# options for star 1
iprofile1 = 0 ! 0=Sink,1=Unif,2=Poly,3=Dens,4=KEPL,5=MESA,6=Pie

# options for star 2
iprofile2 = 0 ! 0=Sink,1=Unif,2=Poly,3=Dens,4=KEPL,5=MESA,6=Pie

Then run phantomsetup again to rewrite the required options::

$ ./phantomsetup sim
run phantomsetup again to rewrite the required options
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

which will give::
::

$ ./phantomsetup sim
...
ERROR: hacc1 not found
ERROR: hacc2 not found
2 error(s) during read of setup file: re-writing...
writing setup options file sim.setup
Edit sim.setup and rerun phantomsetup

Run this again to finish the setup::
run phantomsetup again to complete the setup
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

::
./phantomsetup sim

This creates a sim.in file for the main code. You should edit the tmax and dtmax
This creates a sim.in file for the main code.

edit the .in file as desired
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
You should edit the tmax and dtmax
to give the desired finishing time (tmax) and time between snapshots (dtmax)::

cat sim.in
Expand All @@ -66,20 +97,30 @@ to give the desired finishing time (tmax) and time between snapshots (dtmax)::
f_acc = 0.500 ! particles < f_acc*h_acc accreted without checks


After editing the .in file, proceed to run the simulation::
run the simulation
~~~~~~~~~~~~~~~~~~

::

make
./phantom sim

and have a look at the outputs with splash::
have a look at the outputs with splash
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

::

splash sim_0*


Two polytropes in orbit
~~~~~~~~~~~~~~~~~~~~~~~~~
The default is to setup-and-relax two polytropes and place them in orbit. For this edit
your sim.setup file to give::
------------------------
The default is to setup two polytropes and place them in orbit.

Edit .setup file to specify iprofile=2 for both stars
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

::

# options for star 1
iprofile1 = 2 ! 0=Sink,1=Unif,2=Poly,3=Dens,4=KEPL,5=MESA,6=Pie
Expand All @@ -92,9 +133,30 @@ your sim.setup file to give::
Mstar2 = 1.000 ! mass of star2
Rstar2 = 1.000 ! radius of star2

Then run phantomsetup again to rewrite the required options::
run phantomsetup again to rewrite the required options
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

::

$ ./phantomsetup sim

set the relaxation flag to true
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

::
# relaxation options
relax = T ! relax stars into equilibrium
tol_ekin = 1.000E-07 ! tolerance on ekin/epot to stop relaxation
tol_dens = 1.000 ! % error in density to stop relaxation
maxits = 1000 ! maximum number of relaxation iterations

run phantomsetup again
~~~~~~~~~~~~~~~~~~~~~~~

::

$ ./phantomsetup sim

$ ./phantomsetup sim

This time you should see the automated relax-a-star procedure kick in::

Expand All @@ -104,15 +166,27 @@ This time you should see the automated relax-a-star procedure kick in::
Relaxing star: Iter 2/1000, dens error: 10.97%, R*: 0.915 Ekin/Epot: 4.673E-03
...

As previously, you can then just proceed to run the simulation after editing the sim.in file
::
restart the relaxation if needed
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The relaxation process will continue where it left off if you simply
run phantomsetup again, stopping when the criteria above are met::

$ ./phantomsetup sim


run the simulation after editing the sim.in file
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
After successfully completing the setup process, you can proceed
to run your Very Happy Binary Star Simulation^TM::

./phantom sim.in

Two stars from MESA profiles
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
To use stellar profiles from the MESA code, select iprofile=5 in your .setup file
and enter the name of the ascii data file containing the input profile
-----------------------------

select iprofile=5 in your .setup file and enter name of MESA file
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Enter the name of the ascii data file containing the input profile
(most files produced by MESA just work...)::

# options for star 1
Expand All @@ -131,17 +205,17 @@ and enter the name of the ascii data file containing the input profile
Notice that you do not get to set the particle resolution for the second star,
since the mass of the particles is fixed by the mass and particle number in star 1.

Replacing dense stellar cores with sink particles
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Replace dense stellar cores with sink particles (as necessary)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In the options above you have the option to remove the dense core of the star
which causes small timesteps in the code, and replace it with a softened point
mass. The default option for this is isoftcore1=2 and isinkcore1=1.

For more details, see :doc:`Setting up a softened star <softstar>`


Using SETUP=star and moddump_binary
------------------------------------
Setting up and relaxing binary stars (the old fashioned way)
-------------------------------------------------------------
See :doc:`Setting up stars and tidal disruption events <star>` for the older two-step procedure. The options available are
identical, but with a bit more flexibility and without
having to re-run the relaxation procedure over and over again.
4 changes: 3 additions & 1 deletion scripts/buildbot.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -321,9 +321,11 @@ for setup in $listofsetups; do
colour=$white; # white (default)
ncheck=$((ncheck + 1));

mydebug=''
printf "Checking $setup ($target)... ";
if [[ "$component" == "setup" ]]; then
rm -f $phantomdir/bin/phantomsetup;
mydebug='DEBUG=yes' # compile phantomsetup with DEBUG=yes for setup test
#make clean >& /dev/null;
fi
if [[ "$setup" == "blob" ]]; then
Expand All @@ -332,7 +334,7 @@ for setup in $listofsetups; do
else
mynowarn=$nowarn;
fi
make SETUP=$setup $nolibs $mynowarn $maxp $target 1> $makeout 2> $errorlog; err=$?;
make SETUP=$setup $nolibs $mynowarn $maxp $target $mydebug 1> $makeout 2> $errorlog; err=$?;
#--remove line numbers from error log files
sed -e 's/90(.*)/90/g' -e 's/90:.*:/90/g' $errorlog | grep -v '/tmp' > $errorlog.tmp && mv $errorlog.tmp $errorlog;
if [ $err -eq 0 ]; then
Expand Down
3 changes: 3 additions & 0 deletions src/main/bondiexact.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,9 @@ module bondiexact
real, public :: rcrit = 5.
real, private :: rhocrit = 1.

logical, public :: iswind = .false.
integer, public :: isol = 0

private

contains
Expand Down
1 change: 0 additions & 1 deletion src/main/boundary.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,7 +138,6 @@ subroutine cross_boundary(isperiodic,xyz,ncross)
endif
endif

return
end subroutine cross_boundary

!-----------------------------------------------------------------------
Expand Down
31 changes: 24 additions & 7 deletions src/main/checksetup.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,12 +48,13 @@ subroutine check_setup(nerror,nwarn,restart)
use centreofmass, only:get_centreofmass
use options, only:ieos,icooling,iexternalforce,use_dustfrac,use_hybrid
use io, only:id,master
use externalforces, only:accrete_particles,accradius1,iext_star,iext_corotate
use externalforces, only:accrete_particles,update_externalforce,accradius1,iext_star,iext_corotate
use timestep, only:time
use units, only:G_is_unity,get_G_code
use boundary, only:xmin,xmax,ymin,ymax,zmin,zmax
use boundary_dyn, only:dynamic_bdy,adjust_particles_dynamic_boundary
use nicil, only:n_nden
use metric_tools, only:imetric,imet_minkowski
integer, intent(out) :: nerror,nwarn
logical, intent(in), optional :: restart
integer :: i,nbad,itype,iu,ndead
Expand Down Expand Up @@ -95,11 +96,11 @@ subroutine check_setup(nerror,nwarn,restart)
nwarn = nwarn + 1
endif
#endif
if (hfact < 1. .or. hfact /= hfact) then
if (hfact < 1. .or. isnan(hfact)) then
print*,'ERROR: hfact = ',hfact,', should be >= 1'
nerror = nerror + 1
endif
if (polyk < 0. .or. polyk /= polyk) then
if (polyk < 0. .or. isnan(polyk)) then
print*,'ERROR: polyk = ',polyk,', should be >= 0'
nerror = nerror + 1
endif
Expand Down Expand Up @@ -300,12 +301,18 @@ subroutine check_setup(nerror,nwarn,restart)
!
!--check for particles placed inside accretion boundaries
!
if (iexternalforce > 0 .and. .not.dorestart) then
if (iexternalforce > 0 .and. .not.dorestart .and. (.not.(gr .and. imetric==imet_minkowski))) then
call update_externalforce(iexternalforce,time,0.)
nbad = 0
!$omp parallel do default(none) &
!$omp shared(npart,xyzh,massoftype,time,iexternalforce) &
!$omp private(i,accreted) &
!$omp reduction(+:nbad)
do i=1,npart
call accrete_particles(iexternalforce,xyzh(1,i),xyzh(2,i),xyzh(3,i),xyzh(4,i),massoftype(1),time,accreted)
if (accreted) nbad = nbad + 1
enddo
!$omp end parallel do
if (nbad > 0) then
print*,'WARNING: ',nbad,' of ',npart,' particles setup within the accretion boundary'
nwarn = nwarn + 1
Expand All @@ -314,7 +321,7 @@ subroutine check_setup(nerror,nwarn,restart)
hi = 0.
call accrete_particles(iexternalforce,0.,0.,0.,hi,massoftype(1),time,accreted)
!--if so, check for unresolved accretion radius
if (accreted .and. accradius1 < 0.5*hmin) then
if (accreted .and. accradius1 < 0.5*hmin .and. accradius1 > 0.) then
print*,'WARNING: accretion radius is unresolved by a factor of hmin/racc = ',hmin/accradius1
print*,'(this will cause the code to run needlessly slow)'
nwarn = nwarn + 1
Expand Down Expand Up @@ -450,13 +457,18 @@ subroutine check_NaN(npart,array,label,nerror)
integer :: nbad,i

nbad = 0
!$omp parallel do default(none) schedule(static) &
!$omp shared(npart,array,label) &
!$omp private(i) &
!$omp reduction(+:nbad)
do i=1,npart
!--check for NaNs in xyzh
if (any(isnan(array(:,i)))) then
if (nbad < 10) print*,'NaN in '//trim(label)//' : ', i
nbad = nbad + 1
endif
enddo
!$omp end parallel do
if (nbad > 0) then
print*,'ERROR: NaN in '//trim(label)//' on ',nbad,' of ',npart,' particles'
nerror = nerror + 1
Expand Down Expand Up @@ -844,14 +856,18 @@ subroutine check_for_identical_positions(npart,xyzh,nbad)
!
allocate(index(npart))
call indexxfunc(npart,r2func,xyzh,index)

!
! check for identical positions. Stop checking as soon as non-identical
! positions are found.
!
nbad = 0
itypei = igas
itypej = igas
!$omp parallel do default(none) &
!$omp shared(npart,xyzh,index,maxphase,maxp,iphase) &
!$omp firstprivate(itypei,itypej) &
!$omp private(i,j,dx,dx2) &
!$omp reduction(+:nbad)
do i=1,npart
if (.not.isdead_or_accreted(xyzh(4,index(i)))) then
j = i+1
Expand All @@ -864,7 +880,7 @@ subroutine check_for_identical_positions(npart,xyzh,nbad)
dx2 = dot_product(dx,dx)
if (dx2 < epsilon(dx2) .and. itypei==itypej) then
nbad = nbad + 1
if (nbad <= 100) then
if (nbad <= 10) then
print*,'WARNING: particles of same type at same position: '
print*,' ',index(i),':',xyzh(1:3,index(i))
print*,' ',index(j),':',xyzh(1:3,index(j))
Expand All @@ -874,6 +890,7 @@ subroutine check_for_identical_positions(npart,xyzh,nbad)
enddo
endif
enddo
!$omp end parallel do

deallocate(index)

Expand Down

0 comments on commit 78e4948

Please sign in to comment.