impossibly slow simulations during hyperbolic star collisions

[Francyrad]

Dear users

I wanted to try the new feature to try hyperbolic star collisions like in Starsmasher. However, there are issues: the simulations are impossibly slow, and I hope for my errors.

I attach both the setup and the sim parameters:

# input file for binary setup routines

# units
           mass_unit =      solarm    ! mass unit (e.g. solarm,jupiterm,1e6*solarm)
           dist_unit =      au    ! distance unit (e.g. au,pc,kpc,0.1pc)

# options for star 1
           iprofile1 =           2    ! 0=Sink,1=Unif,2=Poly,3=Dens,4=KEPL,5=MESA,6=Pie
              Mstar1 =       1.000    ! mass of star 1 [Msun]
              Rstar1 =       1.000    ! radius of star1 [Rsun]
                 np1 =        100000    ! number of particles

# options for star 2
           iprofile2 =           2    ! 0=Sink,1=Unif,2=Poly,3=Dens,4=KEPL,5=MESA,6=Pie
              Mstar2 =       0.56    ! mass of star 2 [Msun]
              Rstar2 =       0.50    ! radius of star2 [Rsun]

# orbit settings
                   a =      -0.117    ! semi-major axis (e.g. 1 au) or period (e.g. 10*days)
                 ecc =       1.020    ! eccentricity
                 inc =         90.    ! inclination (deg)
                   O =         90.    ! position angle of ascending node (deg)
                   w =        180.    ! argument of periapsis (deg)
                   f =       -160.    ! initial true anomaly (180=apoastron)
            corotate =           F    ! set stars in corotation

# 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

the sym parameters:

# Runtime options file for Phantom, written 25/02/2024 19:17:22.4
# Options not present assume their default values
# This file is updated automatically after a full dump

# job name
             logfile =  mysim01.log   ! file to which output is directed
            dumpfile =  mysim_00000   ! dump file to start from

# options controlling run time and input/output
                tmax =         10.    ! end time
               dtmax =       1.000    ! time between dumps
                nmax =          -1    ! maximum number of timesteps (0=just get derivs and stop)
                nout =          -1    ! write dumpfile every n dtmax (-ve=ignore)
           nmaxdumps =          -1    ! stop after n full dumps (-ve=ignore)
            twallmax =      000:00    ! maximum wall time (hhh:mm, 000:00=ignore)
           dtwallmax =      024:00    ! maximum wall time between dumps (hhh:mm, 000:00=ignore)
           nfulldump =          10    ! full dump every n dumps
            iverbose =           0    ! verboseness of log (-1=quiet 0=default 1=allsteps 2=debug 5=max)

# options controlling accuracy
              C_cour =       0.300    ! Courant number
             C_force =       0.250    ! dt_force number
                tolv =   1.000E-02    ! tolerance on v iterations in timestepping
               hfact =       1.200    ! h in units of particle spacing [h = hfact(m/rho)^(1/3)]
                tolh =   1.000E-04    ! tolerance on h-rho iterations
       tree_accuracy =       0.500    ! tree opening criterion (0.0-1.0)

# options controlling hydrodynamics, shock capturing
               alpha =       0.000    ! MINIMUM shock viscosity parameter
            alphamax =       1.000    ! MAXIMUM shock viscosity parameter
              alphau =       1.000    ! shock conductivity parameter
                beta =       2.000    ! beta viscosity
        avdecayconst =       0.100    ! decay time constant for viscosity switches

# options controlling damping
               idamp =           0    ! artificial damping of velocities (0=off, 1=constant, 2=star, 3=disc)

# options controlling equation of state
                ieos =           2    ! eqn of state (1=isoth;2=adiab;3=locally iso;8=barotropic)
                  mu =       2.381    ! mean molecular weight
        ipdv_heating =           1    ! heating from PdV work (0=off, 1=on)
      ishock_heating =           1    ! shock heating (0=off, 1=on)

# options controlling cooling
              C_cool =       0.050    ! factor controlling cooling timestep
            icooling =           0    ! cooling function (0=off, 1=cooling library (step), 2=cooling library (force),3=Gammie, 5,6=KI02, 7=powerlaw)

# options controlling sink particles
       icreate_sinks =           0    ! allow automatic sink particle creation
     h_soft_sinksink =       0.000    ! softening length between sink particles
               f_acc =       0.800    ! particles < f_acc*h_acc accreted without checks
      r_merge_uncond =       0.000    ! sinks will unconditionally merge within this separation
        r_merge_cond =       0.000    ! sinks will merge if bound within this radius

# options relating to external forces
      iexternalforce =           0    ! 1=star,2=coro,3=bina,4=prdr,5=toru,6=toys,7=exte,8=spir,9=Lens,10=dens,11=Eins,

# options controlling physical viscosity
           irealvisc =           0    ! physical viscosity type (0=none,1=const,2=Shakura/Sunyaev)
          shearparam =       0.100    ! magnitude of shear viscosity (irealvisc=1) or alpha_SS (irealvisc=2)
            bulkvisc =       0.000    ! magnitude of bulk viscosity

# options for injecting/removing particles
               rkill =      -1.000    ! deactivate particles outside this radius (<0 is off)

# gravitational waves
                  gw =           F    ! calculate gravitational wave strain

The simulation works at full power on Apple Silicon (M1Pro), and it's impossibly slow:

t =  0.45005E-05 dt =  4.529E-06 (courant), np = 156000
 t =  0.90291E-05 dt =  4.557E-06 (courant)
 t =  0.13586E-04 dt =  4.585E-06 (courant)
 t =  0.18170E-04 dt =  4.612E-06 (courant)
 t =  0.22783E-04 dt =  4.640E-06 (courant)
 t =  0.27423E-04 dt =  4.668E-06 (courant)
 t =  0.32091E-04 dt =  4.694E-06 (courant)
 t =  0.36786E-04 dt =  4.721E-06 (courant)
 t =  0.41507E-04 dt =  4.748E-06 (courant)
 t =  0.46255E-04 dt =  4.775E-06 (courant)
 t =  0.51030E-04 dt =  4.801E-06 (courant)
 t =  0.55831E-04 dt =  4.828E-06 (courant)
 t =  0.60659E-04 dt =  4.855E-06 (courant)
 t =  0.65514E-04 dt =  4.882E-06 (courant)
 t =  0.70396E-04 dt =  4.908E-06 (courant)

Am I doing something wrong? I remember once (older version before the hyperbolic implementation) that i ran a binary example with eccentricity = 1 and the simulation was done in a matter of minutes with around 200k particles... Now it can't even run fast 2 identical stars with an eccentricity of 0.

Best

Francesco

[danieljprice]

the limitation of the setup you are using is that you have to set the true anomaly (f, in degrees) to specify the position of the star prior to the encounter. The value you have chosen probably puts the star close to -infinity. Would suggest you have a look at the distance between the two stars from your initial setup and adjust the value of f accordingly

[Francyrad]

The star is not closed to infinity, they are at a distance of 30 solar radius...

The problem is also present when I try to relax the star running the code. The repository that I'm using is some months old, tomorrow I will try the newest and let you know. In case, can you please try to run the code in your laptop? That situation is very weird :(

Thank you for the availability

[Francyrad]

Here i'm again... I installed the new version of the code and the problem is still present: the simulation is impossibly slow using both polytropes and MESA. I have no idea why, i'm sure that months ago my simulation was super fast..

I attach here the setup and the sim setup, am i doing something wrong?

As you can see here, one of the star is not close to infinity:

I have no idea where the problem is, maybe in the setup? I just followed the guide exporting:

SYSTEM = gfortran
export OMP_SCHEDULE="dynamic"
export OMP_STACKSIZE=512M
ulimit -s unlimited

@danieljprice can you please try to run this setup? I have no idea where the problem is:

# input file for binary setup routines

# units
           mass_unit =      solarm    ! mass unit (e.g. solarm,jupiterm,1e6*solarm)
           dist_unit =      au    ! distance unit (e.g. au,pc,kpc,0.1pc)

# options for star 1
           iprofile1 =           2    ! 0=Sink,1=Unif,2=Poly,3=Dens,4=KEPL,5=MESA,6=Pie
              Mstar1 =       1.000    ! mass of star 1 [Msun]
              Rstar1 =       1.000    ! radius of star1 [Rsun]
                 np1 =        100000    ! number of particles

# options for star 2
           iprofile2 =           2    ! 0=Sink,1=Unif,2=Poly,3=Dens,4=KEPL,5=MESA,6=Pie
              Mstar2 =       0.55   ! mass of star 2 [Msun]
              Rstar2 =       0.55    ! radius of star2 [Rsun]

# orbit settings
                   a =      -0.720    ! semi-major axis (e.g. 1 au) or period (e.g. 10*days)
                 ecc =       1.030    ! eccentricity
                 inc =         0.    ! inclination (deg)
                   O =         90.    ! position angle of ascending node (deg)
                   w =        180.    ! argument of periapsis (deg)
                   f =       -100.    ! initial true anomaly (180=apoastron)
            corotate =           F    ! set stars in corotation

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

Thank you for help and availability

Francesco