Add METISSE for stellar evolution

.gitignore

@@ -455,3 +455,5 @@ packages/dist
 
 # standard directory for a virtualenv
 env/
+src/amuse_metisse/data
+src/amuse_metisse/src/METISSE

src/amuse/community/metisse/packages/amuse-metisse.amuse_deps

@@ -0,0 +1 @@
+fortran python cmake install download

src/amuse/community/metisse/packages/amuse-metisse/metisse

@@ -0,0 +1 @@
+../../

src/amuse/community/metisse/packages/amuse-metisse/pyproject.toml

@@ -0,0 +1,44 @@
+[project]
+name = "amuse-metisse"
+dynamic = ["version"]
+requires-python = ">=3.7"
+dependencies = [
+    "amuse-framework"
+]
+
+[build-system]
+requires = ["hatchling", "versioningit"]
+build-backend = "hatchling.build"
+
+[tool.hatch.version]
+source = "versioningit"
+
+[tool.versioningit.next-version]
+method = "smallest"
+
+[tool.versioningit.format]
+distance = "{next_version}.dev{distance}+{vcs}{rev}"
+dirty = "{base_version}+d{build_date:%Y%m%d}"
+distance-dirty = "{next_version}.dev{distance}+{vcs}{rev}.d{build_date:%Y%m%d}"
+
+[tool.hatch.build]
+skip-excluded-dirs = true
+
+[tool.hatch.build.targets.wheel]
+include = ["metisse/**/*.py"]
+exclude = [
+    "metisse/packages",
+    "metisse/support",
+    "metisse/src",
+    "metisse/tests"
+    ]
+artifacts = ["metisse/metisse_worker"]
+
+[tool.hatch.build.targets.wheel.sources]
+"metisse" = "amuse/community/metisse"
+
+[tool.pytest.ini_options]
+pythonpath = ["metisse/tests/", "../../../../../tests"]
+
+testpaths = ["metisse/tests"]
+

src/amuse_metisse/Makefile

@@ -0,0 +1,94 @@
+ifneq (,$(filter-out clean distclean, $(MAKECMDGOALS)))
+include support/config.mk
+endif
+
+
+# Detecting dependencies
+support/config.mk:
+	cd support && ./configure
+
+
+# Downloading the code
+##### Remove if not needed #####
+VERSION = 121257ede5486660d9eb5fbb23b2c377f02ae465
+
+metisse.tar.gz:
+	##### Modify URL as needed #####
+	$(DOWNLOAD) https://github.com/TeamMETISSE/METISSE/archive/$(VERSION).tar.gz >$@
+
+PATCHES := $(file < patches/series)
+PATCHES := $(patsubst %,%,$(PATCHES))    # replace newlines with spaces
+
+src/metisse-$(VERSION): metisse.tar.gz
+	##### Modify as needed #####
+	tar xf $<
+	mv metisse-$(VERSION) src
+	##### Apply patches here, if needed #####
+	#patch <patch.patch
+	#for p in $(PATCHES) ; do patch -p1 <patches/$$p ; done
+
+
+# Building the code into a static library
+DEPFLAGS += $(STOPCOND_CFLAGS) $(STOPCONDMPI_CFLAGS) $(AMUSE_MPI_CFLAGS)
+DEPFLAGS += $(FORSOCKETS_CFLAGS) $(SIMPLE_HASH_CFLAGS) $(G6LIB_CFLAGS)
+DEPFLAGS += $(OPENMP_FFLAGS)
+FCFLAGS += $(DEPFLAGS) -ffpe-summary=all -fbacktrace -fcheck=all -g -Wall
+
+LDLIBS += -lm $(STOPCOND_LIBS) $(STOPCONDMPI_LIBS) $(AMUSE_MPI_LIBS)
+LDLIBS += $(FORSOCKETS_LIBS)
+
+CODELIB := libmetisse.a
+
+##### Modify build command as needed #####
+.PHONY: src/METISSE/$(CODELIB)
+src/METISSE/$(CODELIB): | src/METISSE/make
+	$(MAKE) -C src/METISSE/make libmetisse
+
+
+# Building the workers
+metisse_worker.f90: interface.py
+	amusifier --type=f90 interface.py MetisseInterface -o $@
+
+storage.o: storage.f90
+	$(MPIFC) -c -o $@ $(FCFLAGS) $<
+
+metisse_worker.o: metisse_worker.f90 storage.o
+	$(MPIFC) -c -o $@ $(FCFLAGS) $<
+
+metisse_worker: interface.o metisse_worker.o storage.o src/METISSE/$(CODELIB)
+	$(MPIFC) -o $@ $(LDFLAGS) $^ $(LDLIBS)
+
+%.o: %.f90
+	$(MPIFC) -c -o $@ $(FCFLAGS) -Isrc/METISSE/make $<
+
+
+# Which packages contain which workers?
+amuse-metisse_contains: metisse_worker
+
+# Building and installing packages
+develop-%: %_contains
+	support/shared/uninstall.sh $*
+	python -m pip install -e packages/$*
+
+install-%: %_contains
+	support/shared/uninstall.sh $*
+	python -m pip install packages/$*
+
+package-%: %_contains
+	python3 -m pip install -vv --no-cache-dir --no-deps --no-build-isolation --prefix ${PREFIX} packages/$*
+
+test-%:
+	cd packages/$* && pytest
+
+
+# Cleaning up
+.PHONY: clean
+clean:
+	$(MAKE) -C src clean
+	rm -rf *.o *.mod *worker*
+
+.PHONY: distclean
+distclean: clean
+	rm -f support/config.mk support/config.log support/config.status
+	rm -rf support/autom4te.cache
+

src/amuse_metisse/README

@@ -0,0 +1,13 @@
+
+To update/replace the current version of SSE, just unpack a new
+tarball in the sse directory.  The MUSE module uses the SSE
+distribution unchanged, and relies only on the calling sequences to
+zcnsts(), evolv1(), star(), and deltat().
+
+Scientific papers written using this module should reference the paper
+
+"Comprehensive analytic formulae for stellar evolution as a function 
+ of mass and metallicity" 
+ Hurley J.R., Pols O.R., Tout C.A., 2000, MNRAS, 315, 543 
+
+(See sse/README_SSE for more details.)

src/amuse_metisse/__init__.py

@@ -0,0 +1,5 @@
+"""
+Interface for METISSE
+"""
+from .interface import MetisseInterface
+from .interface import Metisse

src/amuse_metisse/example_metisse.py

@@ -0,0 +1,185 @@
+import logging
+import argparse
+
+import numpy as np
+
+from amuse.units import units, nbody_system
+from amuse.datamodel import Particles
+
+from amuse_metisse import Metisse
+
+from amuse.ic.kroupa import new_kroupa_mass_distribution
+# from amuse.ic.plummer import new_plummer_model
+
+from amuse.io import write_set_to_file
+
+logger = logging.getLogger("amuse")
+# logger.setLevel(logging.DEBUG)
+# logging.basicConfig(level=logging.DEBUG)
+
+
+def new_argument_parser():
+    result = argparse.ArgumentParser(
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+    )
+    result.add_argument("--seed", type=int, default=127)
+    result.add_argument("-n", "--number_of_stars", type=int, default=10000)
+    result.add_argument(
+        "-t",
+        "--time_end",
+        type=units.Myr,
+        default=100.0 | units.Myr,
+        help="The time to evolve to",
+    )
+    result.add_argument(
+        "-s",
+        "--number_of_steps",
+        type=int,
+        default=100,
+        help="Number of time steps (in logspace)",
+    )
+    result.add_argument(
+        "-m",
+        "--mass_min",
+        type=units.MSun,
+        default=0.1 | units.MSun,
+        help="The minimum mass",
+    )
+    result.add_argument(
+        "-M",
+        "--mass_max",
+        type=units.MSun,
+        default=100.0 | units.MSun,
+        help="The maximum mass",
+    )
+    result.add_argument(
+        "--metallicity_dir",
+        type=str,
+        default="/Users/rieder/Code/UvA/Toonen/tres3.0/amuse/src/amuse_metisse/data/Hydrogen",
+        help="The metallicity directory",
+    )
+    result.add_argument(
+        "--metallicity_dir_he",
+        type=str,
+        default="/Users/rieder/Code/UvA/Toonen/tres3.0/amuse/src/amuse_metisse/data/Helium",
+        help="The metallicity-He directory",
+    )
+    return result
+
+
+def setup_metisse(metallicity_dir=None, metallicity_dir_he=None, **kwargs):
+    # instance = Metisse(redirection="file")
+    # instance = Metisse(redirection="none")
+    instance = Metisse()
+
+    instance.parameters.metallicity_dir = metallicity_dir
+    instance.parameters.metallicity_dir_he = metallicity_dir_he
+
+    instance.parameters.wd_mass_scheme = "Modified_mestel"
+    instance.parameters.bhns_mass_scheme = "Belczynski2008"
+    instance.parameters.initial_metallicity = 0.02
+
+    # instance.commit_parameters()
+    return instance
+
+
+def test_metisse_sun(**kwargs):
+    instance = setup_metisse(**kwargs)
+    star = Particles(1)
+    # star.mass = 4.7893208794726441 | units.MSun
+    star.mass = 1.0 | units.MSun
+
+    stars_in_metisse = instance.particles.add_particles(star)
+    print(instance.parameters)
+    assert stars_in_metisse[0].mass == 1.0 | units.MSun
+    # assert stars_in_metisse[0].mass == 4.7893208794726441 | units.MSun
+    print(stars_in_metisse[0])
+    print("Evolving...")
+    instance.evolve_model(1000.0 | units.yr)
+    # instance.evolve_one_step(1)
+    print(stars_in_metisse[0])
+    print("Done")
+    instance.stop()
+
+
+def test_metisse_twostars(**kwargs):
+    instance = setup_metisse(**kwargs)
+    star = Particles(2)
+    star.mass = [0.3, 2.5] | units.MSun
+
+    stars_in_metisse = instance.particles.add_particles(star)
+    print(stars_in_metisse[0])
+    print("Evolving...")
+    instance.evolve_model(1000.0 | units.yr)
+    print(stars_in_metisse.stellar_type)
+    print("Done")
+    instance.stop()
+
+def test_metisse_kroupa(**kwargs):
+    instance = setup_metisse(**kwargs)
+    number_of_stars = 1000
+    star = Particles(number_of_stars)
+    star.mass = new_kroupa_mass_distribution(number_of_stars, mass_min=0.3 | units.MSun, mass_max=10.0 | units.MSun)
+
+    stars_in_metisse = instance.particles.add_particles(star)
+    print(stars_in_metisse[0])
+    print("Evolving...")
+    instance.evolve_model(1000.0 | units.yr)
+    print(stars_in_metisse)
+    print("Done")
+    instance.stop()
+
+
+def evolve_stars_with_metisse(stars, age, **kwargs):
+    instance = setup_metisse(**kwargs)
+    stars_in_metisse = instance.particles.add_particles(stars)
+    instance.evolve_model(age)
+    stars_after_evolution = stars_in_metisse.copy()
+    instance.stop()
+    return stars_after_evolution
+
+
+def evolve_stars_metisse(
+    number_of_stars=100,
+    time_end=100.0 | units.Myr,
+    number_of_steps=100,
+    mass_min=0.1 | units.MSun,
+    mass_max=100.0 | units.MSun,
+    **kwargs
+):
+    mass = np.logspace(
+        np.log10(mass_min.value_in(units.MSun)),
+        np.log10(mass_max.value_in(units.MSun)),
+        number_of_stars
+    ) | units.MSun
+    stars = Particles(number_of_stars)
+    stars.mass = mass
+
+    times = np.logspace(
+        np.log10((time_end / (number_of_steps**1.5)).value_in(units.Myr)),
+        np.log10(time_end.value_in(units.Myr)),
+        number_of_steps
+    ) | units.Myr
+    i = 0
+    for  time in times:
+        print(f"Evolving to time: {time}")
+        stars_metisse = evolve_stars_with_metisse(stars, time, **kwargs)
+        stars_metisse.age = time
+        write_set_to_file(
+            stars_metisse, f"stars_metisse_{i:06d}.amuse"
+        )
+        i += 1
+
+
+def main():
+    args = new_argument_parser().parse_args()
+    np.random.seed(args.seed)
+
+    # test_metisse_sun()
+    # test_metisse_twostars()
+    # test_metisse_kroupa()
+    evolve_stars_metisse(start=0, **vars(args))
+
+
+if __name__ == "__main__":
+    main()