What is the correct way of resuming a simulation with SeBa?

[rieder]

When I stop and later resume a simulation that includes stellar evolution with SeBa, it seems that these stars start again from age 0 when I re-add them to SeBa. How should I rewrite this script so that the stars properly resume their evolution, and that their ages (which aren't equal for all the stars) are properly taken into account?

from amuse.units import units
from amuse.ext.masc import new_star_cluster
from amuse.community.seba.interface import Seba


stars = new_star_cluster(number_of_stars=100)

evo = Seba()
stars_in_seba = evo.particles.add_particles(stars)
code_to_model = stars_in_seba.new_channel_to(stars)
evo.evolve_model(1.0 | units.Myr)
code_to_model.copy()
evo.stop()

print(
    "The average age of stars is %s (± %s)"
    % (stars.age.mean(), stars.age.std())
)

evo = Seba()
stars_in_seba = evo.particles.add_particles(stars)
print(
    "The average age of stars after restarting is %s (± %s)"
    % (stars_in_seba.age.mean(), stars_in_seba.age.std())
)
evo.stop()

[rieder]

The only workaround I can think of is to replay the stellar evolution history completely, which wouldn't be great... I'm hoping there's a better alternative. Any suggestions @spzwart, @silviatoonen?

[silviatoonen]

For now in SeBa, this is the only option. However, if all your stars are single stars (no binaries), it should take little computational efforts to jump to any specific time in the evolution of the stars.

[rieder]

Thanks. For now, I'm dealing with single stars so I can use this workaround. But it would be great if it were possible to resume a calculation - at some point I might want to include binaries :).

For completeness, this is a script that includes the workaround:

import numpy
from amuse.units import units
from amuse.datamodel import Particles
from amuse.ext.masc import new_star_cluster
from amuse.community.seba.interface import Seba


epochs = [0.79, 0.89, 0.94, 0.99] | units.Myr
stars = Particles()
evo = Seba()
offset = None

for time in epochs:
    print("Evolving to time %s" % time)
    if offset is None:
        offset = time
    else:
        evo.evolve_model(time-offset)
    new_stars = new_star_cluster(
        number_of_stars=100,
        lower_mass_limit=2 | units.MSun
    )
    new_stars.birth_mass = new_stars.mass
    new_stars.birth_time = time
    stars.add_particles(new_stars)
    evo.particles.add_particles(new_stars)
    evo.particles.new_channel_to(stars).copy()

evo.stop()

evo = Seba()
stars_with_original_mass = stars.copy()
stars_with_original_mass.mass = stars_with_original_mass.birth_mass
epochs, indices = numpy.unique(
    stars_with_original_mass.birth_time,
    return_inverse=True,
)
offset = None
for i, time in enumerate(epochs):
    if offset is None:
        offset = time
    else:
        evo.evolve_model(time-offset)
    evo.particles.add_particles(
        stars_with_original_mass[indices == i],
    )

[rieder]

For now in SeBa, this is the only option. However, if all your stars are single stars (no binaries), it should take little computational efforts to jump to any specific time in the evolution of the stars.

Unfortunately, at a certain time I will have very many stars, all born at different times. Getting back to the right evolutionary state is more time-consuming than I'd like.

Would it be possible to save the complete state of a star in SeBa, so that it can resume from that state? What would be needed to make this work?

[rieder]

This would be especially of importance when at some point I include binary stars in my simulations, of course.

[silviatoonen]

Hmm, I'd say theoretically if the complete state (whatever that means, something like total mass, relative mass, age, core-mass, stellar-type) is saved, it should be possible. At the moment that infrastructure within SeBa doesn't exist.

[silviatoonen]

After some digging, here's the gist on what is needed: Yes, SeBa can start a simulation with any specific stellar type. This is already possible when using SeBa as a stand alone code. At this moment the implementation in the AMUSE interface is missing -> see line 413 in new_particle() in interface.cc : add_star (.., .., start_type, ..., ...) where start_type is defined on line 18 as static stellar_type start_type = Main_Sequence; What is needed is the option to new_particle to specify the stellar type (note that the stellar types in SeBa are slightly different (fewer) than those of AMUSE).

Besides this one should specify the following parameters: relative_mass envelope_mass core_mass if on the AGB: COcore_mass

current_time relative_age There are different ways of accomplishing this which all need some programming (e.g. adjusting add_star, or making functions which which an amuse user can change the parameters at run time)

[rieder]

Great! I’ll see if I can add the necessary functions for this.

[rieder]

Related question, while I am at it: is it possible to have individual metallicities for the stars? That's currently not allowed by the AMUSE interface but I don't see why SeBa couldn't do it.

[rieder]

(it looks possible so I will include this in my updates)

[silviatoonen]

SeBa-wise it should be fine. Sounds like an Amuse interface thing.

[rieder]

cool. I think it's easiest to add optional parameters to addstar in SeBa, since you'd only want to set these when adding a star. I'll submit a PR to SeBa.

[rieder]

I don't completely understand t_cur and t_current in SeBa's addstar (these appear to share the same value). Currently, when adding a star to SeBa from AMUSE these would always be set to the current time of SeBa (which starts at 0). When restarting, should this be the star's age instead? Or the time the star was born plus the star's age?

[silviatoonen]

For SeBa, ages are always measured from the zero-age main-sequence of that star. aka SeBa doesn't care when the star was born. The thing to pay attention to is the difference between current_time and relative_age, the former being the actual time passed since the zero-age main-sequence, the latter being a measure for how far a long a stellar evolution track the star is.

[stale[bot]]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed in 28 days if no further activity occurs. Thank you for your contributions.