diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 8098ff3743c52eb01cb9537a59ab734e07929c58..fbaeba529ff1208399d417e7e00419057632cc9b 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -35,6 +35,5 @@ before_script:
test:
script:
- - cd test
- - ./getCoolingTable.sh
- - ./test_cooling.py
\ No newline at end of file
+ - cd tests
+ - ./run_tests.sh
\ No newline at end of file
diff --git a/docs/RST/examples.rst b/docs/RST/examples.rst
index 0912efdbcafde647c6059e69b21fbedb34a6d76f..ef5fc94927e18969ac0b05d00edefa2c20675d1b 100644
--- a/docs/RST/examples.rst
+++ b/docs/RST/examples.rst
@@ -8,3 +8,4 @@ If you wish to generate them, you will need to first install PySWIFTsim and then
:maxdepth: 2
examples/cooling.rst
+ examples/stellar.rst
diff --git a/docs/RST/examples/cooling.rst b/docs/RST/examples/cooling.rst
index 1cc6f4db5c44493cdf7792754084760d5b89954b..f72b3f8287868a6c5780fe8f44fca7620ddcffd7 100644
--- a/docs/RST/examples/cooling.rst
+++ b/docs/RST/examples/cooling.rst
@@ -6,21 +6,25 @@ Cooling Rate
~~~~~~~~~~~~
This example is generated in ``examples/cooling_rate`` with the Grackle cooling.
-The plot are the same than in [Smith2016]_.
+The plot are the same than in [Smith2016]_ and uses a solar metallicity.
The code has two different modes 1D or 2D.
In the first mode, the code generates particles at different temperatures but same density and computes the cooling rate associated to theses conditions.
-.. image:: https://obswww.unige.ch/~lhausamm/PySWIFTsim/examples/cooling_rate.png
+.. image:: https://obswww.unige.ch/lastro/projects/Clastro/PySWIFTsim/examples/cooling_rate.png
In the second mode, the code generates a grid of particles at different density and temperatures and then computes the cooling.
+I suppose that in [Smith2016]_, they use the non equilibrium mode of Grackle to generate this graph.
+It requires a careful initialization of the element fractions which is not easy due to the non monotonic behavior of the cooling function.
+For simplicity, only Grackle with the equilibrium mode is shown in this documentation.
-.. image:: https://obswww.unige.ch/~lhausamm/PySWIFTsim/examples/cooling_rate_2d.png
+
+.. image:: https://obswww.unige.ch/lastro/projects/Clastro/PySWIFTsim/examples/cooling_rate_2d.png
Cooling Box
~~~~~~~~~~~
-.. image:: https://obswww.unige.ch/~lhausamm/PySWIFTsim/examples/cooling_box.png
+.. image:: https://obswww.unige.ch/lastro/projects/Clastro/PySWIFTsim/examples/cooling_box.png
diff --git a/docs/RST/examples/stellar.rst b/docs/RST/examples/stellar.rst
new file mode 100644
index 0000000000000000000000000000000000000000..94e1ac8ee917b03ae0ee1fd50d3886c07f312707
--- /dev/null
+++ b/docs/RST/examples/stellar.rst
@@ -0,0 +1,25 @@
+Stellar Examples
+================
+
+
+Initial Mass Function
+~~~~~~~~~~~~~~~~~~~~~
+
+This example is generated in ``examples/initial_mass_function`` with GEAR's stellar model.
+The initial mass function used in GEAR follows [Kroupa2001]_.
+
+.. image:: https://obswww.unige.ch/lastro/projects/Clastro/PySWIFTsim/examples/initial_mass_function.png
+
+
+Lifetime
+~~~~~~~~
+
+This example is generated in ``examples/lifetime`` with GEAR's stellar model.
+The lifetime in GEAR is a quadratic in metallicity and mass (see equation 3.32 in [Poirier2004]_)
+
+
+.. image:: https://obswww.unige.ch/lastro/projects/Clastro/PySWIFTsim/examples/lifetime.png
+
+
+.. [Kroupa2001] `On the variation of the Initial Mass Function <https://arxiv.org/abs/astro-ph/0009005>`_
+.. [Poirier2004] `Étude de l'évolution chimique et dynamique d'objets proto-galactiques: Application à l'évolution des galaxies spirales <https://obswww.unige.ch/~revaz/PyChem/_downloads/ff03ba21204a3f36fa4ec885abd545ae/ThesePoirier.pdf>`_
diff --git a/docs/RST/libstellar.rst b/docs/RST/libstellar.rst
new file mode 100644
index 0000000000000000000000000000000000000000..1a4950817b25d46c1aa53ab9bb12f8b8365669e8
--- /dev/null
+++ b/docs/RST/libstellar.rst
@@ -0,0 +1,11 @@
+Libstellar
+==========
+
+.. automodule:: pyswiftsim.libstellar
+ :members:
+ :undoc-members:
+ :private-members:
+ :special-members:
+ :show-inheritance:
+ :inherited-members:
+
diff --git a/docs/RST/pyswiftsim.rst b/docs/RST/pyswiftsim.rst
index 0e77bb58ed4a0765bdd419362edebdda15f60178..d4b10ad0820d4c377a87ea1537b938d325e8b92e 100644
--- a/docs/RST/pyswiftsim.rst
+++ b/docs/RST/pyswiftsim.rst
@@ -5,7 +5,4 @@ PySWIFTsim
:members:
:undoc-members:
:private-members:
- :special-members:
- :show-inheritance:
- :inherited-members:
diff --git a/docs/index.rst b/docs/index.rst
index d11523e49bcf892c03137a516830300e753817cc..a799ccce26e3c2226423ef2d296779ad3c3bdaf2 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -14,6 +14,7 @@ The aim is to provide a python tool that benefits from HPC code and allows to di
RST/pyswiftsim.rst
RST/libcooling.rst
+ RST/libstellar.rst
RST/examples.rst
diff --git a/examples/cooling_box/cooling.yml b/examples/cooling_box/cooling.yml
deleted file mode 100644
index 415f22ce7fadc74f042ad91d30079aa8ce0ffbea..0000000000000000000000000000000000000000
--- a/examples/cooling_box/cooling.yml
+++ /dev/null
@@ -1,28 +0,0 @@
-# Define the system of units to use internally.
-InternalUnitSystem:
- UnitMass_in_cgs: 1.989e43 # Grams
- UnitLength_in_cgs: 3.085e21 # Centimeters
- UnitVelocity_in_cgs: 1e5 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
-
-# Parameters for the hydrodynamics scheme
-SPH:
- resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
-
-# Cooling with Grackle 2.0
-GrackleCooling:
- CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table
- WithUVbackground: 1 # Enable or not the UV background
- Redshift: 0 # Redshift to use (-1 means time based redshift)
- WithMetalCooling: 1 # Enable or not the metal cooling
- ProvideVolumetricHeatingRates: 0 # User provide volumetric heating rates
- ProvideSpecificHeatingRates: 0 # User provide specific heating rates
- SelfShieldingMethod: 0 # Grackle (<= 3) or Gear self shielding method
- OutputMode: 0 # Write in output corresponding primordial chemistry mode
- MaxSteps: 1000
- ConvergenceLimit: 1e-2
-
-GearChemistry:
- InitialMetallicity: 0.01295
\ No newline at end of file
diff --git a/examples/cooling_box/plot_cooling.py b/examples/cooling_box/plot_cooling.py
index 13c2e4e514fd14b4f749a8a7d7f39a1e19414088..18bbb26903a25e87667887dc5beeacc8f3feb976 100644
--- a/examples/cooling_box/plot_cooling.py
+++ b/examples/cooling_box/plot_cooling.py
@@ -18,7 +18,7 @@
# ########################################################################
from pyswiftsim.libcooling import coolingRate
-from pyswiftsim import downloadCoolingTable
+import pyswiftsim
import numpy as np
import matplotlib.pyplot as plt
@@ -28,24 +28,28 @@ from time import strftime, gmtime
plt.rc('text', usetex=True)
-downloadCoolingTable()
+pyswiftsim.downloadCoolingTable()
# PARAMETERS
# swift params filename
filename = "cooling.yml"
+# use cosmology?
with_cosmo = False
+# include metal cooling?
+with_metals = 1
+
# particle data
# in atom/cc
-part_density = np.array([1.])
-part_temperature = np.array([1e4]) # in K
+part_density = np.array([0.1])
+part_temperature = np.array([1e6]) # in K
gamma = 5./3.
# time in Myr
-t_end = 1.
-N_time = 100000
+t_end = 1e2
+N_time = 1000
# SCRIPT
@@ -86,51 +90,68 @@ def plot_solution(energy, data, us):
unit_length = float(us["UnitLength_in_cgs"])
unit_vel = float(us["UnitVelocity_in_cgs"])
unit_mass = float(us["UnitMass_in_cgs"])
+ unit_temp = float(us["UnitTemp_in_cgs"])
unit_time = unit_length / unit_vel
# change units => cgs
rho *= unit_mass / unit_length**3
- energy *= unit_length**2 / unit_time**2
+ m_p = constants.m_p.to("g") / (unit_mass * units.g)
+ rho = part_density * unit_length**3 * m_p
+ d["density"] = rho.to("")
+
+ unit_energy = unit_mass * unit_length**2
+ unit_energy /= unit_time**2
+ k_B = constants.k_B.to("erg / K") / (unit_energy / unit_temp)
+ k_B *= units.K / units.erg
+ T = unit_temp * energy / k_B
+ T *= (gamma - 1.) * m_p
- time *= unit_time
+ myr = units.s / units.Myr
+ time *= unit_time * float(myr.to(""))
# do plot
plt.figure()
date = strftime("%d %b %Y", gmtime())
plt.title("Date: {}".format(date))
- plt.plot(time, energy)
- plt.xlabel("Time [s]")
- plt.ylabel("Energy [erg]")
- plt.show()
+ plt.plot(time, T)
+ plt.xlabel("Time [Myr]")
+ plt.ylabel("Temperature [K]")
+ ax = plt.gca()
+ ax.set_xscale("log")
+ ax.set_yscale("log")
+ ax.set_xlim([1e-2, 1e2])
-def getParser(filename):
- with open(filename, "r") as stream:
- stream = yaml.load(stream)
- return stream
+ plt.show()
if __name__ == "__main__":
- us = getParser(filename)["InternalUnitSystem"]
+ overwrite = {
+ "GrackleCooling": {
+ "WithMetalCooling": with_metals,
+ }
+ }
+ with pyswiftsim.ParameterFile(overwrite) as filename:
+ us = pyswiftsim.parseYamlFile(filename)["InternalUnitSystem"]
- d = generate_initial_condition(us)
+ d = generate_initial_condition(us)
- t = d["time"]
- dt = t[1] - t[0]
+ t = d["time"]
+ dt = t[1] - t[0]
- # du / dt
- print("Computing cooling...")
+ # du / dt
+ print("Computing cooling...")
- N = len(t)
- energy = np.zeros(t.shape)
- energy[0] = d["energy"]
- rho = np.array(d["density"], dtype=np.float32)
- for i in range(N-1):
- ene = np.array([energy[i]], dtype=np.float32)
- rate = coolingRate(filename, rho, ene, with_cosmo, dt)
- energy[i+1] = energy[i] + rate * dt
+ N = len(t)
+ energy = np.zeros(t.shape)
+ energy[0] = d["energy"]
+ rho = np.array(d["density"], dtype=np.float32)
+ for i in range(N-1):
+ ene = np.array([energy[i]], dtype=np.float32)
+ rate = coolingRate(filename, rho, ene, with_cosmo, dt)
+ energy[i+1] = energy[i] + rate * dt
- plot_solution(energy, d, us)
+ plot_solution(energy, d, us)
diff --git a/examples/cooling_rate/cooling.yml b/examples/cooling_rate/cooling.yml
deleted file mode 100644
index 166fa0401a5870b941e503152d9843a8e2417e69..0000000000000000000000000000000000000000
--- a/examples/cooling_rate/cooling.yml
+++ /dev/null
@@ -1,41 +0,0 @@
-# Define the system of units to use internally.
-InternalUnitSystem:
- UnitMass_in_cgs: 1.989e43 # Grams
- UnitLength_in_cgs: 3.085e21 # Centimeters
- UnitVelocity_in_cgs: 1e5 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
-
-# Parameters for the hydrodynamics scheme
-SPH:
- resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
- minimal_temperature: 0 # (Optional) Minimal temperature (in internal units) allowed for the gas particles. Value is ignored if set to 0.
-
-# Cosmological parameters
-Cosmology:
- h: 0.6777 # Reduced Hubble constant
- a_begin: 0.0078125 # Initial scale-factor of the simulation
- a_end: 1.0 # Final scale factor of the simulation
- Omega_m: 0.307 # Matter density parameter
- Omega_lambda: 0.693 # Dark-energy density parameter
- Omega_b: 0.0455 # Baryon density parameter
- Omega_r: 0. # (Optional) Radiation density parameter
- w_0: -1.0 # (Optional) Dark-energy equation-of-state parameter at z=0.
- w_a: 0. # (Optional) Dark-energy equation-of-state time evolution parameter.
-
-# Cooling with Grackle 2.0
-GrackleCooling:
- CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table
- WithUVbackground: 1 # Enable or not the UV background
- Redshift: 0 # Redshift to use (-1 means time based redshift)
- WithMetalCooling: 1 # Enable or not the metal cooling
- ProvideVolumetricHeatingRates: 0 # User provide volumetric heating rates
- ProvideSpecificHeatingRates: 0 # User provide specific heating rates
- SelfShieldingMethod: 0 # Grackle (<= 3) or Gear self shielding method
- ConvergenceLimit: 1e-4
- MaxSteps: 20000
- Omega: 0.8
-
-GearChemistry:
- InitialMetallicity: 0.01295
\ No newline at end of file
diff --git a/examples/cooling_rate/plot_cooling.py b/examples/cooling_rate/plot_cooling.py
index 55d21ee1790177aa1755a7cc2815cdb26414502a..884f0592da44acc0444e14e41ea63dcba67bf4eb 100644
--- a/examples/cooling_rate/plot_cooling.py
+++ b/examples/cooling_rate/plot_cooling.py
@@ -18,31 +18,34 @@
# ########################################################################
from pyswiftsim.libcooling import coolingRate
-from pyswiftsim import downloadCoolingTable
+import pyswiftsim
from copy import deepcopy
import numpy as np
import matplotlib.pyplot as plt
from astropy import units, constants
-import yaml
from time import strftime, gmtime
+import sys
plt.rc('text', usetex=True)
-downloadCoolingTable()
+pyswiftsim.downloadCoolingTable()
# PARAMETERS
+
+colormap = "RdBu"
+
# cosmology
with_cosmo = False
-# swift params filename
-filename = "cooling.yml"
+# include metals?
+with_metals = 1
# hydrogen mass fraction
h_mass_frac = 0.76
# density in atom / cm3
-N_rho = 100
+N_rho = int(sys.argv[-1])
if N_rho == 1:
density = np.array([1.])
else:
@@ -110,12 +113,6 @@ def generate_initial_condition(us):
return d
-def getParser(filename):
- with open(filename, "r") as stream:
- stream = yaml.load(stream)
- return stream
-
-
def plot_solution(rate, data, us):
print("Plotting solution")
energy = data["energy"]
@@ -169,7 +166,8 @@ def plot_solution(rate, data, us):
N_levels = 100
levels = np.linspace(_min, _max, N_levels)
plt.figure()
- plt.contourf(rho, T, cooling_length, levels)
+ plt.contourf(rho, T, cooling_length, levels,
+ cmap=plt.get_cmap(colormap))
plt.xlabel("Density [atom/cm3]")
plt.ylabel("Temperature [K]")
@@ -190,17 +188,24 @@ def plot_solution(rate, data, us):
if __name__ == "__main__":
- parser = getParser(filename)
- us = parser["InternalUnitSystem"]
+ overwrite = {
+ "GrackleCooling": {
+ "WithMetalCooling": with_metals,
+ }
+ }
+ with pyswiftsim.ParameterFile(overwrite) as filename:
+
+ parser = pyswiftsim.parseYamlFile(filename)
+ us = parser["InternalUnitSystem"]
- d = generate_initial_condition(us)
+ d = generate_initial_condition(us)
- # du / dt
- print("Computing cooling...")
+ # du / dt
+ print("Computing cooling...")
- rate = coolingRate(filename,
- d["density"].astype(np.float32),
- d["energy"].astype(np.float32),
- with_cosmo, d["time_step"])
+ rate = coolingRate(filename,
+ d["density"].astype(np.float32),
+ d["energy"].astype(np.float32),
+ with_cosmo, d["time_step"])
- plot_solution(rate, d, us)
+ plot_solution(rate, d, us)
diff --git a/examples/cooling_rate/run.sh b/examples/cooling_rate/run.sh
index ca4ed090d7161d18ce520af16a81caaaebd04cf3..ce3d510430474c3986b463878c3edff2ba02d212 100644
--- a/examples/cooling_rate/run.sh
+++ b/examples/cooling_rate/run.sh
@@ -1,3 +1,4 @@
#!/bin/bash
-./plot_cooling.py
+./plot_cooling.py 1
+./plot_cooling.py 100
diff --git a/examples/initial_mass_function/plot_initial_mass_function.py b/examples/initial_mass_function/plot_initial_mass_function.py
new file mode 100644
index 0000000000000000000000000000000000000000..bc653edc70a36756b7f65f206fd4e5e02cd56129
--- /dev/null
+++ b/examples/initial_mass_function/plot_initial_mass_function.py
@@ -0,0 +1,53 @@
+#!/usr/bin/env python3
+# ########################################################################
+# This file is part of PYSWIFT.
+# Copyright (c) 2018 Loic Hausammann (loic.hausammann@epfl.ch)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+# ########################################################################
+
+from pyswiftsim import libstellar
+import pyswiftsim
+
+import numpy as np
+import matplotlib.pyplot as plt
+plt.rc('text', usetex=True)
+
+N = 100
+
+
+if __name__ == "__main__":
+
+ with pyswiftsim.ParameterFile() as filename:
+
+ parser = pyswiftsim.parseYamlFile(filename)
+ imf_parser = parser["GEARInitialMassFunction"]
+
+ mass_limits = list(imf_parser["mass_limits"])
+ mass_min = np.log10(mass_limits[0])
+ mass_max = np.log10(mass_limits[-1])
+
+ # du / dt
+ print("Computing initial mass function...")
+
+ mass = np.logspace(mass_min, mass_max, N, dtype=np.float32)
+ imf = libstellar.initialMassFunction(filename, mass)
+
+ # plot IMF
+ plt.loglog(mass, imf)
+
+ plt.xlabel("Mass [Msun]")
+ plt.ylabel("Mass fraction")
+
+ plt.show()
diff --git a/examples/initial_mass_function/run.sh b/examples/initial_mass_function/run.sh
new file mode 100644
index 0000000000000000000000000000000000000000..78b3761d0e26f631942cf03a5def8cf55104f51d
--- /dev/null
+++ b/examples/initial_mass_function/run.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+./plot_initial_mass_function.py
diff --git a/examples/lifetime/plot_lifetime.py b/examples/lifetime/plot_lifetime.py
new file mode 100644
index 0000000000000000000000000000000000000000..7696311e3cfd048d2dd000bebbbfc4389d448125
--- /dev/null
+++ b/examples/lifetime/plot_lifetime.py
@@ -0,0 +1,63 @@
+#!/usr/bin/env python3
+# ########################################################################
+# This file is part of PYSWIFT.
+# Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+# ########################################################################
+
+from pyswiftsim import libstellar
+import pyswiftsim
+
+import numpy as np
+import matplotlib.pyplot as plt
+plt.rc('text', usetex=True)
+
+
+N = 1000
+solar_metallicity = 0.02041
+
+
+if __name__ == "__main__":
+ with pyswiftsim.ParameterFile() as filename:
+
+ parser = pyswiftsim.parseYamlFile(filename)
+ # Get masses
+ imf = parser["GEARInitialMassFunction"]
+ mass_limits = imf["mass_limits"]
+ mass_min = np.log10(mass_limits[0])
+ mass_max = np.log10(mass_limits[-1])
+
+ mass = np.logspace(mass_min, mass_max, N, dtype=np.float32)
+
+ # get metallicity
+ Z = solar_metallicity * np.ones(N, dtype=np.float32)
+
+ print("Computing lifetime...")
+
+ lifetime = libstellar.lifetimeFromMass(filename, mass, Z)
+ mass_from_lifetime = libstellar.massFromLifetime(filename, lifetime, Z)
+
+ lifetime /= 1e6
+
+ plt.loglog(mass, lifetime, label="Lifetime from mass")
+ plt.loglog(mass_from_lifetime, lifetime, "--",
+ label="Mass from lifetime")
+
+ plt.xlabel("Mass [Msun]")
+ plt.ylabel("Lifetime [Myr]")
+
+ plt.legend()
+
+ plt.show()
diff --git a/examples/lifetime/run.sh b/examples/lifetime/run.sh
new file mode 100644
index 0000000000000000000000000000000000000000..a180bae803170bc4e66db1432a2b367dffcdfe36
--- /dev/null
+++ b/examples/lifetime/run.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+./plot_lifetime.py
diff --git a/pyswiftsim/__init__.py b/pyswiftsim/__init__.py
index 8c6a5748ae78c94a83aee7c1d808247ac8fd102e..7852bc713c33d950db8883ca686b71f12e998b4f 100644
--- a/pyswiftsim/__init__.py
+++ b/pyswiftsim/__init__.py
@@ -18,9 +18,15 @@
import os
import urllib.request
+from tempfile import NamedTemporaryFile
+import yaml
def downloadCoolingTable():
+ """
+ Download the cooling table if it does not already exist in the current
+ repository.
+ """
url = "http://virgodb.cosma.dur.ac.uk/"
url += "swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5"
filename = "CloudyData_UVB=HM2012.h5"
@@ -28,3 +34,283 @@ def downloadCoolingTable():
if not os.path.isfile(filename):
# Download the file from `url` and save it locally under `file_name`:
urllib.request.urlretrieve(url, filename)
+
+
+def downloadYieldsTable():
+ """
+ Download the yields table if it does not already exist in the current
+ repository.
+ """
+ url = "https://obswww.unige.ch/"
+ url += "lastro/projects/Clastro/PySWIFTsim/"
+ filename = "chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5"
+
+ if not os.path.isfile(filename):
+ # Download the file from `url` and save it locally under `file_name`:
+ urllib.request.urlretrieve(url, filename)
+
+
+def parseYamlFile(filename):
+ """
+ Parse a yaml file and return it as a dictionaryw.
+
+ Parameters
+ ==========
+
+ filename: str
+ The filename of the yaml file
+
+ Returns
+ =======
+
+ yaml: dict
+ The parameters
+
+ Examples
+ ========
+
+ >>> yaml = parseYamlFile("swift.yml")
+ >>> us = yaml["InternalUnitSystem"]
+
+ Everything in the dictionary is a string, therefore you will need to
+ convert them manually.
+
+ >>> unit_mass = float(us["UnitMass_in_cgs"])
+ """
+ with open(filename, "r") as stream:
+ stream = yaml.load(stream)
+ return stream
+
+
+class ParameterFile:
+ """
+ This structure generates a temporary parameter file that is automatically
+ removed. If the code crashes before the destruction of the parameter file,
+ the parameter file will not be removed, but it should be taken care
+ of by the operating system.
+
+ Parameters
+ ==========
+
+ overwrite: dict
+ Allow to overwrite some default parameters or to add new parameters.
+
+ existing: str
+ Filename of an existing parameter file.
+
+ Examples
+ ========
+
+ >>> overwrite = {"InternalUnitSystem": {"UnitMass_in_cgs": 1.}}
+ >>> with ParameterFile(overwrite=overwrite) as filename:
+ >>> params = parseYamlFile(filename)
+ >>> unit_mass = params["InternalUnitSystem"]["UnitMass_in_cgs"]
+ >>> print(unit_mass)
+ 1.
+
+
+ When leaving the ``with`` environment, the file is deleted.
+ """
+ default_parameters = {
+ "InternalUnitSystem": {
+ "UnitMass_in_cgs": 1.989e43, # 1e10 Msun
+ "UnitLength_in_cgs": 3.085e21, # Kpc
+ "UnitVelocity_in_cgs": 1e5, # km / s
+ "UnitCurrent_in_cgs": 1, # Amperes
+ "UnitTemp_in_cgs": 1, # Kelvin
+ },
+
+ "SPH": {
+ "resolution_eta": 1.2348,
+ "CFL_condition": 0.1,
+ "minimal_temperature": 0
+ },
+
+ "Cosmology": {
+ "h": 0.6777,
+ "a_begin": 0.0078125,
+ "a_end": 1.0,
+ "Omega_m": 0.307,
+ "Omega_lambda": 0.693,
+ "Omega_b": 0.0455,
+ "Omega_r": 0.,
+ "w_0": -1.0,
+ "w_a": 0.
+ },
+
+ "GrackleCooling": {
+ "CloudyTable": "CloudyData_UVB=HM2012.h5",
+ "WithUVbackground": 1,
+ "Redshift": 0,
+ "WithMetalCooling": 1,
+ "ProvideVolumetricHeatingRates": 0,
+ "ProvideSpecificHeatingRates": 0,
+ "SelfShieldingMethod": 0,
+ "ConvergenceLimit": 1e-2,
+ "MaxSteps": 20000,
+ "Omega": 0.8
+ },
+
+ "GearChemistry": {
+ "InitialMetallicity": 0.01295,
+ },
+
+ "GEARInitialMassFunction": {
+ "number_function_part": 4,
+ "exponents": [0.7, -0.8, -1.7, -1.3],
+ "mass_limits": [0.05, 0.08, 0.5, 1.0, 50.]
+ },
+
+ "GEARFeedback": {
+ "SupernovaeEnergy_erg": 1e51,
+ "ThermalTime_Myr": 5,
+ "ChemistryTable": "chemistry.hdf5",
+ },
+
+ "GEARLifetime": {
+ "quadratic": [-40.110, 5.509, 0.7824],
+ "linear": [141.929, -15.889, -3.2557],
+ "constant": [-261.365, 17.073, 9.8661],
+ }
+ }
+ """
+ Dictionnary containing all the default parameters to write in the
+ parameter file.
+ """
+
+ def _write(self, f, txt):
+ """
+ Encode and then write the string into a stream.
+
+ Parameters
+ ==========
+
+ f: stream (e.g. an opened file)
+ The stream where to write the information
+
+ txt: str
+ The text to write in the steam
+
+ """
+ txt = txt.encode()
+ f.write(txt)
+
+ def _writeParameters(self, f, d, overwrite):
+ """
+ Writes the parameters file.
+
+ Parameters
+ ==========
+
+ f: stream (e.g. opened file)
+ The stream where to write the parameter file.
+
+ d: dict
+ The dictionary containing the default parameters to write.
+
+ overwrite: dict
+ The dictionary containing the parameters to overwrite.
+ """
+ _format = " {}: {}\n"
+ for sec in d.keys():
+ p = d[sec]
+ # write section
+ txt = "{}:\n".format(sec)
+ self._write(f, txt)
+ for param in p.keys():
+ value = p[param]
+
+ # test if overwritten
+ test = overwrite is not None and sec in overwrite
+ test = test and param in overwrite[sec]
+ if test:
+ value = overwrite[sec][param]
+ # write the parameter
+ txt = _format.format(param, value)
+ self._write(f, txt)
+
+ # check if in the overwrite dictionary
+ if overwrite is None or sec not in overwrite:
+ continue
+
+ # ensure that all the elements in overwrite are written
+ for param in overwrite[sec].keys():
+ # check if already done
+ if param in p:
+ continue
+
+ # write the parameter
+ txt = _format.format(param, overwrite[sec][param])
+ self._write(f, txt)
+
+ # add a space between sections
+ print()
+
+ if overwrite is None:
+ return
+
+ # now do all the sections that do not exist in default
+ for sec in overwrite.keys():
+ # check if already done
+ if sec in d:
+ continue
+
+ # write section
+ txt = "{}:\n".format(sec)
+ self._write(f, txt)
+ for param in overwrite[sec].keys():
+ # write parameter
+ txt = _format.format(param, overwrite[sec][param])
+ self._write(f, txt)
+
+ # add space between sections
+ print()
+
+ def __init__(self, overwrite=None, existing=None):
+ """
+ This function Initializes the structure, but does not create the file.
+
+ Parameters
+ ==========
+
+ overwrite: dict
+ Allow to overwrite some default parameters or to add new parameters.
+
+ existing: str
+ Filename of an existing parameter file.
+
+ """
+ if overwrite is not None and existing is not None:
+ raise Exception("Cannot overwrite default parameters and"
+ " use existing parameter file")
+ self.overwrite = overwrite
+ self.d = self.default_parameters
+ self.existing = existing
+
+ def __enter__(self):
+ """
+ Create the parameter file.
+
+ Returns
+ =======
+
+ filename: str
+ The filename of the parameter file.
+ """
+ if self.existing is not None:
+ return self.existing
+
+ f = NamedTemporaryFile(delete=False)
+ self.filename = f.name
+
+ self._writeParameters(f, self.default_parameters, self.overwrite)
+
+ f.close()
+ return self.filename
+
+ def __exit__(self, *args):
+ """
+ Remove the temporary file.
+ """
+ if self.existing is None:
+ os.remove(self.filename)
diff --git a/scripts/pyswift_cooling_rate b/scripts/pyswift_cooling_rate
new file mode 100644
index 0000000000000000000000000000000000000000..262629c6debc640b90821d828f9190ef5b67d872
--- /dev/null
+++ b/scripts/pyswift_cooling_rate
@@ -0,0 +1,278 @@
+#!/usr/bin/env python3
+# ########################################################################
+# This file is part of PYSWIFT.
+# Copyright (c) 2018 Loic Hausammann (loic.hausammann@epfl.ch)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+# ########################################################################
+
+from pyswiftsim.libcooling import coolingRate
+import pyswiftsim
+
+from copy import deepcopy
+import numpy as np
+import matplotlib.pyplot as plt
+from astropy import units, constants
+import yaml
+from time import strftime, gmtime
+from argparse import ArgumentParser
+import os
+
+plt.rc('text', usetex=True)
+
+# PARAMETERS
+
+colormap = "RdBu"
+
+# cosmology
+with_cosmo = False
+
+
+def parseArguments():
+ parser = ArgumentParser(description="Plot the cooling rate")
+
+ parser.add_argument(
+ "--download_table", dest="table",
+ action="store_true",
+ help="Download the default cooling table")
+
+ parser.add_argument(
+ "--params", dest="params",
+ type=str, default=None,
+ help="SWIFT parameters file")
+
+ parser.add_argument(
+ "--T_min", dest="Tmin",
+ type=float, default=10.,
+ help="Minimal temperature to plot")
+
+ parser.add_argument(
+ "--T_max", dest="Tmax",
+ type=float, default=1e9,
+ help="Maximal temperature to plot")
+
+ parser.add_argument(
+ "--N_T", dest="N_T",
+ type=int, default=1000,
+ help="Number of points for the temperature")
+
+ parser.add_argument(
+ "--rho_min", dest="rho_min",
+ type=float, default=1e-6,
+ help="Minimal density to plot")
+
+ parser.add_argument(
+ "--rho_max", dest="rho_max",
+ type=float, default=1e4,
+ help="Maximal density to plot (only used if more than 1 point)")
+
+ parser.add_argument(
+ "--N_rho", dest="N_rho",
+ type=int, default=1,
+ help="Number of points for the density")
+
+ parser.add_argument(
+ "--gamma", dest="gamma",
+ type=float, default=5. / 3.,
+ help="Adiabatic index")
+
+ parser.add_argument(
+ "--h_frac", dest="h_mass_frac",
+ type=float, default=0.76,
+ help="hydogen mass fraction")
+
+ args = parser.parse_args()
+
+ return args
+
+
+def mu(T, h_mass_frac):
+ # define constants
+ T_transition = 1e4
+ mu_neutral = 4. / (1. + 3. * h_mass_frac)
+ mu_ion = 4. / (8. - 5. * (1 - h_mass_frac))
+
+ # Find correct mu
+ ind = T > T_transition
+ mu = np.zeros(T.shape)
+ mu[ind] = mu_ion
+ mu[~ind] = mu_neutral
+ return mu
+
+
+# SCRIPT
+def generate_initial_condition(us, opt):
+ print("Generating default initial conditions")
+ # get units
+ unit_length = float(us["UnitLength_in_cgs"])
+ unit_vel = float(us["UnitVelocity_in_cgs"])
+ unit_mass = float(us["UnitMass_in_cgs"])
+ unit_temp = float(us["UnitTemp_in_cgs"])
+
+ # density in atom / cm3
+ N_rho = opt.N_rho
+ if N_rho == 1:
+ density = np.array([opt.rho_min])
+ else:
+ rho_min = np.log10(opt.rho_min)
+ rho_max = np.log10(opt.rho_max)
+ density = np.logspace(rho_min, rho_max, N_rho)
+
+ # temperature in K
+ N_T = opt.N_T
+ T_min = np.log10(opt.Tmin)
+ T_max = np.log10(opt.Tmax)
+ temperature = np.logspace(T_min, T_max, N_T)
+
+ # time step in s
+ dt = units.Myr * 1e-8
+ dt = dt.to("s") / units.s
+
+ d = {}
+ # generate grid
+ rho, T = np.meshgrid(density, temperature)
+ rho = deepcopy(rho.reshape(-1))
+ T = T.reshape(-1)
+ d["temperature"] = T
+
+ # Deal with units
+ m_p = constants.m_p.to("g") / (unit_mass * units.g)
+ rho = rho * unit_length**3 * m_p
+ d["density"] = rho.to("")
+
+ unit_time = unit_length / unit_vel
+ unit_energy = unit_mass * unit_length**2
+ unit_energy /= unit_time**2
+ k_B = constants.k_B.to("erg / K") / (unit_energy / unit_temp)
+ k_B *= units.K / units.erg
+ energy = k_B * T / (unit_temp * mu(T, opt.h_mass_frac))
+ energy /= (opt.gamma - 1.) * m_p
+ d["energy"] = energy.to("")
+
+ time_step = dt / unit_time
+ d["time_step"] = time_step
+
+ return d
+
+
+def getParser(filename):
+ with open(filename, "r") as stream:
+ stream = yaml.load(stream)
+ return stream
+
+
+def plot_solution(rate, data, us, opt):
+ print("Plotting solution")
+ energy = data["energy"]
+ rho = data["density"]
+ T = data["temperature"]
+
+ T *= float(us["UnitTemp_in_cgs"])
+
+ # change units => cgs
+ u_m = float(us["UnitMass_in_cgs"])
+ u_l = float(us["UnitLength_in_cgs"])
+ u_v = float(us["UnitVelocity_in_cgs"])
+ u_t = u_l / u_v
+
+ rho *= u_m / u_l**3
+
+ energy *= u_v**2
+
+ rate *= u_v**2 / u_t
+
+ # lambda cooling
+ lam = rate * rho
+
+ N_rho = opt.N_rho
+ N_T = opt.N_T
+
+ # do plot
+ if N_rho == 1:
+ # plot Lambda vs T
+ plt.figure()
+ plt.loglog(T, np.abs(lam),
+ label="SWIFT")
+ plt.xlabel("Temperature [K]")
+ plt.ylabel("$\\Lambda$ [erg s$^{-1}$ cm$^{3}$]")
+
+ else:
+ m_p = constants.m_p.to("g") / units.g
+ rho /= m_p
+
+ shape = [N_rho, N_T]
+ cooling_time = energy / rate
+ cooling_length = np.sqrt(opt.gamma * (opt.gamma-1.) * energy)
+ cooling_length *= cooling_time
+
+ cooling_length = np.log10(np.abs(cooling_length) / units.kpc.to('cm'))
+
+ # reshape
+ rho = rho.reshape(shape)
+ T = T.reshape(shape)
+ energy = energy.reshape(shape)
+ cooling_length = cooling_length.reshape(shape)
+
+ _min = -7
+ _max = 7
+ N_levels = 100
+ levels = np.linspace(_min, _max, N_levels)
+ plt.figure()
+ plt.contourf(rho, T, cooling_length, levels,
+ cmap=plt.get_cmap(colormap))
+ plt.xlabel("Density [atom/cm3]")
+ plt.ylabel("Temperature [K]")
+
+ ax = plt.gca()
+ ax.set_xscale("log")
+ ax.set_yscale("log")
+
+ cbar = plt.colorbar()
+ tc = np.arange(_min, _max, 1.5)
+ cbar.set_ticks(tc)
+ cbar.set_ticklabels(tc)
+ cbar.ax.set_ylabel("Log( Cooling Length / kpc )")
+
+ date = strftime("%d %b %Y", gmtime())
+ plt.title("Date: {}".format(date))
+ plt.show()
+
+
+if __name__ == "__main__":
+
+ opt = parseArguments()
+
+ if opt.table:
+ pyswiftsim.downloadCoolingTable()
+
+ if opt.params is not None and not os.path.isfile(opt.params):
+ raise Exception(
+ "The parameter file '{}' does not exist".format(opt.params))
+
+ with pyswiftsim.ParameterFile(existing=opt.params) as filename:
+
+ parser = getParser(filename)
+ us = parser["InternalUnitSystem"]
+
+ d = generate_initial_condition(us, opt)
+
+ # du / dt
+ print("Computing cooling...")
+
+ rate = coolingRate(filename,
+ d["density"].astype(np.float32),
+ d["energy"].astype(np.float32),
+ with_cosmo, d["time_step"])
+
+ plot_solution(rate, d, us, opt)
diff --git a/scripts/pyswift_initial_mass_function b/scripts/pyswift_initial_mass_function
new file mode 100644
index 0000000000000000000000000000000000000000..803e520495fc9d453963975de2e4100d20374bc5
--- /dev/null
+++ b/scripts/pyswift_initial_mass_function
@@ -0,0 +1,101 @@
+#!/usr/bin/env python3
+# ########################################################################
+# This file is part of PYSWIFT.
+# Copyright (c) 2018 Loic Hausammann (loic.hausammann@epfl.ch)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+# ########################################################################
+
+from pyswiftsim import libstellar
+import pyswiftsim
+
+import numpy as np
+import matplotlib.pyplot as plt
+from argparse import ArgumentParser
+import os
+
+plt.rc('text', usetex=True)
+
+
+def parseArguments():
+ parser = ArgumentParser(description="Plot the initial mass function")
+
+ parser.add_argument(
+ "--download_table", dest="table",
+ action="store_true",
+ help="Download the default chemistry table")
+
+ parser.add_argument(
+ "--params", dest="params",
+ type=str, default=None,
+ help="SWIFT parameters file")
+
+ parser.add_argument(
+ "--N", dest="N",
+ type=int, default=1000,
+ help="Number of points")
+
+ parser.add_argument(
+ "--mass_min", dest="mass_min",
+ type=float, default=None,
+ help="Minimal mass to plot")
+
+ parser.add_argument(
+ "--mass_max", dest="mass_max",
+ type=float, default=None,
+ help="Maximal mass to plot")
+
+ args = parser.parse_args()
+
+ return args
+
+
+if __name__ == "__main__":
+ opt = parseArguments()
+
+ if opt.table:
+ pyswiftsim.downloadChemistryTable()
+
+ if opt.params is not None and not os.path.isfile(opt.params):
+ raise Exception(
+ "The parameter file '{}' does not exist".format(opt.params))
+
+ with pyswiftsim.ParameterFile(existing=opt.params) as filename:
+
+ parser = pyswiftsim.parseYamlFile(filename)
+ imf_parser = parser["GEARInitialMassFunction"]
+
+ mass_limits = list(imf_parser["mass_limits"])
+ if opt.mass_min is not None:
+ mass_min = np.log10(opt.mass_min)
+ else:
+ mass_min = np.log10(mass_limits[0])
+
+ if opt.mass_max is not None:
+ mass_max = np.log10(opt.mass_max)
+ else:
+ mass_max = np.log10(mass_limits[-1])
+
+ print("Computing initial mass function...")
+
+ mass = np.logspace(mass_min, mass_max, opt.N, dtype=np.float32)
+ imf = libstellar.initialMassFunction(filename, mass)
+
+ # plot IMF
+ plt.loglog(mass, imf)
+
+ plt.xlabel("Mass [Msun]")
+ plt.ylabel("Mass fraction")
+
+ plt.show()
diff --git a/scripts/pyswift_lifetime b/scripts/pyswift_lifetime
new file mode 100644
index 0000000000000000000000000000000000000000..0b14c53f589d9a26c21841b4393536a000dc1eec
--- /dev/null
+++ b/scripts/pyswift_lifetime
@@ -0,0 +1,109 @@
+#!/usr/bin/env python3
+# ########################################################################
+# This file is part of PYSWIFT.
+# Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+# ########################################################################
+
+from pyswiftsim import libstellar
+import pyswiftsim
+
+import numpy as np
+import matplotlib.pyplot as plt
+from argparse import ArgumentParser
+import os
+plt.rc('text', usetex=True)
+
+
+def parseArguments():
+ parser = ArgumentParser(description="Plot the initial mass function")
+
+ parser.add_argument(
+ "--download_table", dest="table",
+ action="store_true",
+ help="Download the default chemistry table")
+
+ parser.add_argument(
+ "--params", dest="params",
+ type=str, default=None,
+ help="SWIFT parameters file")
+
+ parser.add_argument(
+ "--N", dest="N",
+ type=int, default=1000,
+ help="Number of points")
+
+ parser.add_argument(
+ "--mass_min", dest="mass_min",
+ type=float, default=None,
+ help="Minimal mass to plot")
+
+ parser.add_argument(
+ "--mass_max", dest="mass_max",
+ type=float, default=None,
+ help="Maximal mass to plot")
+
+ parser.add_argument(
+ "--metallicity", dest="metallicity",
+ type=float, default=0.02041,
+ help="Metallicity of the stars")
+
+ args = parser.parse_args()
+
+ return args
+
+
+if __name__ == "__main__":
+ opt = parseArguments()
+
+ if opt.table:
+ pyswiftsim.downloadChemistryTable()
+
+ if opt.params is not None and not os.path.isfile(opt.params):
+ raise Exception(
+ "The parameter file '{}' does not exist".format(opt.params))
+
+ with pyswiftsim.ParameterFile(existing=opt.params) as filename:
+
+ parser = pyswiftsim.parseYamlFile(filename)
+ imf_parser = parser["GEARInitialMassFunction"]
+
+ mass_limits = list(imf_parser["mass_limits"])
+ if opt.mass_min is not None:
+ mass_min = np.log10(opt.mass_min)
+ else:
+ mass_min = np.log10(mass_limits[0])
+
+ if opt.mass_max is not None:
+ mass_max = np.log10(opt.mass_max)
+ else:
+ mass_max = np.log10(mass_limits[-1])
+
+ mass = np.logspace(mass_min, mass_max, opt.N, dtype=np.float32)
+
+ # get metallicity
+ Z = opt.metallicity * np.ones(opt.N, dtype=np.float32)
+
+ print("Computing lifetime...")
+
+ lifetime = libstellar.lifetimeFromMass(filename, mass, Z)
+ lifetime /= 1e6
+
+ plt.loglog(mass, lifetime)
+
+ plt.xlabel("Mass [Msun]")
+ plt.ylabel("Lifetime [Myr]")
+
+ plt.show()
diff --git a/setup.py b/setup.py
index 1c76bc94e8463de8a2bf4af34c8a9c61a2505365..62a38d6006d0bd59d22d2e1efdd3254612664d7c 100644
--- a/setup.py
+++ b/setup.py
@@ -21,10 +21,6 @@ cflags = [
# disables some warnings due to python
"-Wno-unused-parameter",
"-Wno-maybe-uninitialized",
- "-Wno-strict-prototypes",
- "-Wno-unused-function",
- "-Wno-incompatible-pointer-types",
- "-Wno-missing-field-initializers",
"-fopenmp",
]
@@ -118,15 +114,22 @@ data_files = []
##############
c_src = glob("src/*.c")
-ext_modules = Extension("pyswiftsim.libcooling",
- c_src,
- include_dirs=include,
- libraries=lib,
- library_dirs=lib_dir,
- extra_compile_args=cflags,
- extra_link_args=lflags)
-
-ext_modules = [ext_modules]
+ext_modules = []
+ext_modules.append(Extension("pyswiftsim.libcooling",
+ c_src,
+ include_dirs=include,
+ libraries=lib,
+ library_dirs=lib_dir,
+ extra_compile_args=cflags,
+ extra_link_args=lflags))
+
+ext_modules.append(Extension("pyswiftsim.libstellar",
+ c_src,
+ include_dirs=include,
+ libraries=lib,
+ library_dirs=lib_dir,
+ extra_compile_args=cflags,
+ extra_link_args=lflags))
##############
# data
@@ -138,7 +141,7 @@ data_files = []
# scripts
##############
-list_scripts = []
+list_scripts = glob("scripts/*")
##############
# Setup
diff --git a/src/cooling_wrapper.c b/src/cooling_wrapper.c
index 47dd1ccb52857d1794719eaac9ba2bdcedded1c2..4fa9731fb44d1a66db652c95c553ab590ddc8dd9 100644
--- a/src/cooling_wrapper.c
+++ b/src/cooling_wrapper.c
@@ -66,7 +66,7 @@ void pycooling_set_fractions(struct xpart *xp, PyArrayObject* frac, const int id
* @brief Compute the cooling rate
*
* args is expecting:
- * - a parameter filename,
+ * - the name of the parameter file
* - a numpy array of the densities,
* - a numpy array of the specific energies,
* - an optional int if the cosmology should be used,
@@ -117,44 +117,14 @@ PyObject* pycooling_rate(PyObject* self, PyObject* args) {
size_t N = PyArray_DIM(energy, 0);
- /* Initialize parser */
- struct swift_params params;
- parser_read_file(filename, ¶ms);
+ PYSWIFTSIM_INIT_STRUCTS();
- /* Init unit system */
- struct unit_system us;
- units_init_from_params(&us, ¶ms, "InternalUnitSystem");
-
- /* Init physical constant */
- struct phys_const pconst;
- phys_const_init(&us, ¶ms, &pconst);
-
- /* init cooling */
- if (initialized == 0) {
- // struct cooling_function_data cooling;
- cooling_init_backend(¶ms, &us, &pconst, &cooling);
- initialized = 1;
- }
- else if (initialized == 1) {
- message("WARNING: not reinitializing the cooling, need to be fixed");
- initialized = 2;
- }
+ PYSWIFTSIM_INIT_COOLING();
/* Init chemistry */
struct chemistry_global_data chemistry;
chemistry_init_backend(¶ms, &us, &pconst, &chemistry);
- /* Init cosmo */
- struct cosmology cosmo;
-
- if (with_cosmo)
- cosmology_init(¶ms, &us, &pconst, &cosmo);
- else
- cosmology_init_no_cosmo(&cosmo);
-
- /* Init hydro prop */
- struct hydro_props hydro_props;
- hydro_props_init(&hydro_props, &pconst, &us, ¶ms);
/* Init entropy floor */
struct entropy_floor_properties floor_props;
@@ -177,12 +147,17 @@ PyObject* pycooling_rate(PyObject* self, PyObject* args) {
float u = *(float*) PyArray_GETPTR1(energy, i);
hydro_set_physical_internal_energy(&p, &xp, &cosmo, u);
+ chemistry_first_init_part(&pconst, &us, &cosmo, &chemistry, &p, &xp);
+
if (fractions != NULL)
pycooling_set_fractions(&xp, fractions, i);
- else
- cooling_first_init_part(&pconst, &us, &cosmo, &cooling, &p, &xp);
+ else {
+ /* Need to set the mass in order to estimate the element fractions */
+ hydro_set_mass(&p, p.rho);
+ p.h = 1;
- chemistry_first_init_part(&pconst, &us, &cosmo, &chemistry, &p, &xp);
+ cooling_first_init_part(&pconst, &us, &cosmo, &cooling, &p, &xp);
+ }
hydro_set_physical_internal_energy_dt(&p, &cosmo, 0);
diff --git a/src/cooling_wrapper.h b/src/cooling_wrapper.h
index bd9017f8abb4c31084db4ace60fa52372f943006..e0b97a9c9f70d9033dfaac54ae75bb061f8dfb35 100644
--- a/src/cooling_wrapper.h
+++ b/src/cooling_wrapper.h
@@ -28,5 +28,20 @@ extern int initialized;
PyObject* pycooling_rate(PyObject* self, PyObject* args);
+#define PYSWIFTSIM_INIT_COOLING() \
+ \
+ /* init cooling */ \
+ if (initialized == 0) { \
+ /* struct cooling_function_data cooling; */ \
+ cooling_init_backend(¶ms, &us, &pconst, &cooling); \
+ initialized = 1; \
+ } \
+ else if (initialized == 1) { \
+ message("WARNING: not reinitializing the cooling, need to be fixed"); \
+ initialized = 2; \
+ } \
+ \
+
+
#endif // __PYSWIFTSIM_COOLING_H__
diff --git a/src/libstellar.c b/src/libstellar.c
new file mode 100644
index 0000000000000000000000000000000000000000..fcd1cb304d5207fd409eb9461337e0bac03ac60e
--- /dev/null
+++ b/src/libstellar.c
@@ -0,0 +1,130 @@
+/*******************************************************************************
+ * This file is part of PYSWIFTSIM.
+ * Copyright (c) 2018 loic hausammann (loic.hausammann@epfl.ch)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#include "stellar_evolution_wrapper.h"
+
+#include <Python.h>
+#include <math.h>
+#include <numpy/arrayobject.h>
+
+/* definition of the method table */
+
+static PyMethodDef libstellar_methods[] = {
+
+ {"initialMassFunction", pyinitial_mass_function, METH_VARARGS,
+ "Compute the initial mass function.\n\n"
+ "Parameters\n"
+ "----------\n\n"
+ "filename: str\n"
+ "\t Parameter file\n"
+ "mass: array\n"
+ "\t Mass of the stars in internal units.\n"
+ "\n"
+ "Returns\n"
+ "-------\n"
+ "imf: array\n"
+ "\t The initial mass function.\n"
+ "\n"
+ "Examples\n"
+ "--------\n"
+ ">>> mass = np.array([1, 10])\n"
+ ">>> filename = 'params.yml'\n"
+ ">>> imf = libstellar.initialMassFunction(filename, mass)\n"
+
+ },
+
+ {"lifetimeFromMass", pylifetime_from_mass, METH_VARARGS,
+ "Compute the lifetime of a star from its mass and metallicity.\n\n"
+ "Parameters\n"
+ "----------\n\n"
+ "filename: str\n"
+ "\t Parameter file\n"
+ "mass: array\n"
+ "\t Mass of the stars in internal units.\n"
+ "metallicity: array\n"
+ "\t Metallicity of the stars.\n"
+ "\n"
+ "Returns\n"
+ "-------\n"
+ "tau: array\n"
+ "\t The lifetime of the stars.\n"
+ "\n"
+ "Examples\n"
+ "--------\n"
+ ">>> mass = np.array([1, 10])\n"
+ ">>> Z = np.array([0.02, 0.01])\n"
+ ">>> filename = 'params.yml'\n"
+ ">>> tau = libstellar.lifetimeFromMass(filename, mass, Z)\n"
+ },
+
+ {"massFromLifetime", pymass_from_lifetime, METH_VARARGS,
+ "Compute the mass of a star from its lifetime and metallicity.\n\n"
+ "Parameters\n"
+ "----------\n\n"
+ "filename: str\n"
+ "\t Parameter file\n"
+ "lifetime: array\n"
+ "\t Lifetime of the stars in internal units.\n"
+ "metallicity: array\n"
+ "\t Metallicity of the stars.\n"
+ "\n"
+ "Returns\n"
+ "-------\n"
+ "mass: array\n"
+ "\t The mass of the stars.\n"
+ "\n"
+ "Examples\n"
+ "--------\n"
+ ">>> lifetime = np.array([1e7, 1e9])\n"
+ ">>> Z = np.array([0.02, 0.01])\n"
+ ">>> filename = 'params.yml'\n"
+ ">>> mass = libstellar.massFromLifetime(filename, lifetime, Z)\n"
+ },
+
+
+ {NULL, NULL, 0, NULL} /* Sentinel */
+};
+
+
+
+static struct PyModuleDef libstellar_cmodule = {
+ PyModuleDef_HEAD_INIT,
+ "libstellar",
+ "Wrapper around the stellar physics of SWIFT",
+ -1,
+ libstellar_methods,
+ NULL, /* m_slots */
+ NULL, /* m_traverse */
+ NULL, /* m_clear */
+ NULL, /* m_free */
+};
+
+
+PyMODINIT_FUNC PyInit_libstellar(void)
+{
+ PyObject *m;
+
+ /* set time for swift */
+ clocks_set_cpufreq(0);
+
+ import_array();
+ Py_Initialize();
+ m = PyModule_Create(&libstellar_cmodule);
+
+ return m;
+}
diff --git a/src/pyswiftsim_tools.h b/src/pyswiftsim_tools.h
index 567314f028f17bf34e395ea77292763d1373447e..6b75075c21bc10bd6399cce584e128c83adb54ee 100644
--- a/src/pyswiftsim_tools.h
+++ b/src/pyswiftsim_tools.h
@@ -44,4 +44,32 @@ enum error_code {
int pytools_check_array(PyArrayObject *obj, int dim, int type, const char* name);
+#define PYSWIFTSIM_INIT_STRUCTS() \
+ \
+ /* Initialize parser */ \
+ struct swift_params params; \
+ parser_read_file(filename, ¶ms); \
+ \
+ /* Init unit system */ \
+ struct unit_system us; \
+ units_init_from_params(&us, ¶ms, "InternalUnitSystem"); \
+ \
+ /* Init physical constant */ \
+ struct phys_const pconst; \
+ phys_const_init(&us, ¶ms, &pconst); \
+ \
+ /* Init cosmo */ \
+ struct cosmology cosmo; \
+ \
+ if (with_cosmo) \
+ cosmology_init(¶ms, &us, &pconst, &cosmo); \
+ else \
+ cosmology_init_no_cosmo(&cosmo); \
+ \
+ /* Init hydro prop */ \
+ struct hydro_props hydro_props; \
+ hydro_props_init(&hydro_props, &pconst, &us, ¶ms); \
+
+
+
#endif // __PYSWIFTSIM_TOOLS_H__
diff --git a/src/stellar_evolution_wrapper.c b/src/stellar_evolution_wrapper.c
new file mode 100644
index 0000000000000000000000000000000000000000..a990b753d2c05ddc620f2115bce366a75fc4ce09
--- /dev/null
+++ b/src/stellar_evolution_wrapper.c
@@ -0,0 +1,192 @@
+/*******************************************************************************
+ * This file is part of PYSWIFTSIM.
+ * Copyright (c) 2018 loic hausammann (loic.hausammann@epfl.ch)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#include "stellar_evolution_wrapper.h"
+#include "pyswiftsim_tools.h"
+
+/**
+ * @brief Compute the initial mass function.
+ *
+ * args is expecting:
+ * - the name of the parameter file
+ * - a numpy array containing the masses.
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return The initial mass function
+ */
+PyObject* pyinitial_mass_function(PyObject* self, PyObject* args) {
+ import_array();
+
+
+ char *filename;
+
+ PyArrayObject *mass;
+
+ /* parse argument */
+ if (!PyArg_ParseTuple(
+ args, "sO", &filename, &mass))
+ return NULL;
+
+ /* check numpy array */
+ if (pytools_check_array(mass, 1, NPY_FLOAT, "Mass") != SWIFT_SUCCESS)
+ return NULL;
+
+ size_t N = PyArray_DIM(mass, 0);
+
+ int with_cosmo = 0;
+ PYSWIFTSIM_INIT_STRUCTS();
+
+ /* Init stellar physics */
+ struct feedback_props fp;
+ feedback_props_init(&fp, &pconst, &us, ¶ms, &hydro_props, &cosmo);
+
+ /* return object */
+ PyArrayObject *imf = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(mass), PyArray_DIMS(mass), NPY_FLOAT);
+
+ for(size_t i = 0; i < N; i++) {
+ float m = *(float*) PyArray_GETPTR1(mass, i);
+ float *imf_i = (float*) PyArray_GETPTR1(imf, i);
+
+
+ *imf_i = stellar_evolution_get_imf(&fp.stellar_model.imf, m);
+ }
+
+ return (PyObject *) imf;
+
+}
+
+/**
+ * @brief Compute the lifetime of a star.
+ *
+ * args is expecting:
+ * TODO
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return cooling rate
+ */
+PyObject* pylifetime_from_mass(PyObject* self, PyObject* args) {
+ import_array();
+
+ char *filename;
+
+ PyArrayObject *mass;
+ PyArrayObject *metallicity;
+
+ /* parse argument */
+ if (!PyArg_ParseTuple(
+ args, "sOO", &filename, &mass, &metallicity))
+ return NULL;
+
+ /* check numpy array */
+ if (pytools_check_array(mass, 1, NPY_FLOAT, "Mass") != SWIFT_SUCCESS)
+ return NULL;
+ if (pytools_check_array(metallicity, 1, NPY_FLOAT, "Metallicity") != SWIFT_SUCCESS)
+ return NULL;
+
+ if (PyArray_DIM(metallicity, 0) != PyArray_DIM(mass, 0))
+ error("Mass and metallicity should have the same dimension");
+
+
+ size_t N = PyArray_DIM(mass, 0);
+
+ int with_cosmo = 0;
+
+ PYSWIFTSIM_INIT_STRUCTS();
+
+ /* Init stellar physics */
+ struct feedback_props fp;
+ feedback_props_init(&fp, &pconst, &us, ¶ms, &hydro_props, &cosmo);
+
+ /* return object */
+ PyArrayObject *lifetime = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(mass), PyArray_DIMS(mass), NPY_FLOAT);
+
+ for(size_t i = 0; i < N; i++) {
+ float m = *(float*) PyArray_GETPTR1(mass, i);
+ float z = *(float*) PyArray_GETPTR1(metallicity, i);
+ float *tau = (float*) PyArray_GETPTR1(lifetime, i);
+
+
+ *tau = stellar_evolution_get_log_lifetime_from_mass(&fp.stellar_model.lifetime, log10(m), z);
+ *tau = pow(10, *tau);
+ }
+
+ return (PyObject *) lifetime;
+
+}
+
+/**
+ * @brief Compute the mass of a star with a given lifetime.
+ *
+ * args is expecting:
+ * TODO
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return cooling rate
+ */
+PyObject* pymass_from_lifetime(PyObject* self, PyObject* args) {
+ import_array();
+
+ char *filename;
+
+ PyArrayObject *time;
+ PyArrayObject *metallicity;
+
+ /* parse argument */
+ if (!PyArg_ParseTuple(
+ args, "sOO", &filename, &time, &metallicity))
+ return NULL;
+
+ /* check numpy array */
+ if (pytools_check_array(time, 1, NPY_FLOAT, "Lifetime") != SWIFT_SUCCESS)
+ return NULL;
+ if (pytools_check_array(metallicity, 1, NPY_FLOAT, "Metallicity") != SWIFT_SUCCESS)
+ return NULL;
+
+ if (PyArray_DIM(metallicity, 0) != PyArray_DIM(time, 0))
+ error("Lifetime and metallicity should have the same dimension");
+
+
+ size_t N = PyArray_DIM(time, 0);
+
+ int with_cosmo = 0;
+
+ PYSWIFTSIM_INIT_STRUCTS();
+
+ /* Init stellar physics */
+ struct feedback_props fp;
+ feedback_props_init(&fp, &pconst, &us, ¶ms, &hydro_props, &cosmo);
+
+ /* return object */
+ PyArrayObject *mass = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(time), PyArray_DIMS(time), NPY_FLOAT);
+
+ for(size_t i = 0; i < N; i++) {
+ float tau = *(float*) PyArray_GETPTR1(time, i);
+ float z = *(float*) PyArray_GETPTR1(metallicity, i);
+ float *m = (float*) PyArray_GETPTR1(mass, i);
+
+
+ *m = stellar_evolution_get_log_mass_from_lifetime(&fp.stellar_model.lifetime, log10(tau), z);
+ *m = pow(10, *m);
+ }
+
+ return (PyObject *) mass;
+
+}
diff --git a/src/stellar_evolution_wrapper.h b/src/stellar_evolution_wrapper.h
new file mode 100644
index 0000000000000000000000000000000000000000..a76b04d39b1733a8667cc9435859cfeb9fbfbeea
--- /dev/null
+++ b/src/stellar_evolution_wrapper.h
@@ -0,0 +1,29 @@
+/*******************************************************************************
+ * This file is part of PYSWIFTSIM.
+ * Copyright (c) 2018 loic hausammann (loic.hausammann@epfl.ch)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef __PYSWIFTSIM_STELLAR_H__
+#define __PYSWIFTSIM_STELLAR_H__
+
+#include "pyswiftsim_includes.h"
+
+PyObject* pyinitial_mass_function(PyObject* self, PyObject* args);
+PyObject* pylifetime_from_mass(PyObject* self, PyObject* args);
+PyObject* pymass_from_lifetime(PyObject* self, PyObject* args);
+
+#endif // __PYSWIFTSIM_STELLAR_H__
+
diff --git a/test/test_cooling.py b/test/test_cooling.py
deleted file mode 100644
index 94e3a46917f74659c8b8502f2922d316f14d5c41..0000000000000000000000000000000000000000
--- a/test/test_cooling.py
+++ /dev/null
@@ -1,86 +0,0 @@
-#!/usr/bin/env python3
-# ########################################################################
-# This file is part of PYSWIFT.
-# Copyright (c) 2018 Loic Hausammann (loic.hausammann@epfl.ch)
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published
-# by the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-# ########################################################################
-
-from pyswiftsim.libcooling import coolingRate
-from pyswiftsim import downloadCoolingTable
-import numpy as np
-from astropy import units, constants
-import yaml
-
-# define parameters
-filename = "test_cooling.yml"
-# adiabatic index
-gamma = 5. / 3.
-
-# download cooling table
-downloadCoolingTable()
-
-# Read parameter file
-f = open(filename, "r")
-data = yaml.load(f)
-f.close()
-
-# Get the units
-swift_units = data["InternalUnitSystem"]
-# length unit
-u_l = float(swift_units["UnitLength_in_cgs"])
-# velocity unit
-u_v = float(swift_units["UnitVelocity_in_cgs"])
-# mass unit
-u_m = float(swift_units["UnitMass_in_cgs"])
-# time unit
-u_t = u_l / u_v
-# energy unit
-u_e = u_m * u_l**2 / u_t**2
-
-# boltzman constant
-k_B = constants.k_B.to("erg / K") * units.K / (units.erg * u_e)
-k_B = float(k_B)
-# proton mass
-m_p = constants.m_p.to("g") / (units.g * u_m)
-m_p = float(m_p)
-
-
-def internalEnergy(T):
- u = (k_B / m_p) * T
- u /= (gamma - 1.)
- return u
-
-
-# generates data
-# time step in Myr
-dt = 1e-3
-dt = float(units.Myr.to('s')) / u_t
-
-# density in atom / cm3
-density = np.array([1.], dtype=np.float32)
-density *= m_p * u_l**3
-
-# Temperature in K
-T = np.array([1e4], dtype=np.float32)
-u = internalEnergy(T)
-
-# compute the cooling rate
-with_cosmo = False
-rate = coolingRate(filename, density, u, with_cosmo, dt)
-
-# go back to cgs
-rate /= u_e / u_t
-
-print("Cooling done: {} erg/s".format(rate))
diff --git a/test/test_cooling.yml b/test/test_cooling.yml
deleted file mode 100644
index a30e9edf8a5fa99564e288e5029c038af1fb2b32..0000000000000000000000000000000000000000
--- a/test/test_cooling.yml
+++ /dev/null
@@ -1,25 +0,0 @@
-# Define the system of units to use internally.
-InternalUnitSystem:
- UnitMass_in_cgs: 1.989e43 # Grams
- UnitLength_in_cgs: 3.085e21 # Centimeters
- UnitVelocity_in_cgs: 1e5 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
-
-# Cooling with Grackle 3.0
-GrackleCooling:
- CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
- WithUVbackground: 1 # Enable or not the UV background
- Redshift: 0 # Redshift to use (-1 means time based redshift)
- WithMetalCooling: 1 # Enable or not the metal cooling
- ProvideVolumetricHeatingRates: 0 # (optional) User provide volumetric heating rates
- ProvideSpecificHeatingRates: 0 # (optional) User provide specific heating rates
- SelfShieldingMethod: 0 # (optional) Grackle (<= 3) or Gear self shielding method
- OutputMode: 0 # (optional) Write in output corresponding primordial chemistry mode
- MaxSteps: 10000 # (optional) Max number of step when computing the initial composition
- ConvergenceLimit: 1e-2 # (optional) Convergence threshold (relative) for initial composition
-
-# Parameters for the hydrodynamics scheme
-SPH:
- resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
diff --git a/tests/generate_parameter_file.py b/tests/generate_parameter_file.py
new file mode 100644
index 0000000000000000000000000000000000000000..576823abc6a284f8c79281e412f0d0de69dce657
--- /dev/null
+++ b/tests/generate_parameter_file.py
@@ -0,0 +1,85 @@
+#!/usr/bin/env python3
+# ########################################################################
+# This file is part of PYSWIFT.
+# Copyright (c) 2018 Loic Hausammann (loic.hausammann@epfl.ch)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+# ########################################################################
+
+from pyswiftsim import ParameterFile, parseYamlFile
+import os
+
+# defines a few parameter that will be overwritten
+overwrite = {
+ "InternalUnitSystem": {
+ "UnitMass_in_cgs": 1,
+ },
+
+ "test": {
+ "test": -1
+ }
+}
+
+print("Testing the default file")
+
+# Generate and read the default file
+with ParameterFile() as params:
+ default = parseYamlFile(params)
+
+# check if the temporary file is removed
+if os.path.isfile(params):
+ raise Exception("File {} not destroyed".format(params))
+
+
+print("Testing to overwrite the default file")
+
+# Now do the same but with the overwrite feature
+with ParameterFile(overwrite) as params:
+ over = parseYamlFile(params)
+
+# check if the temporary file is removed
+if os.path.isfile(params):
+ raise Exception("File {} not destroyed".format(params))
+
+print("Comparing results")
+
+# Now check that the correct values are found
+for sec in over.keys():
+ params = over[sec]
+
+ for p in params.keys():
+ # First case not in default parameters
+ if sec not in default:
+ assert overwrite[sec][p] == params[p]
+ continue
+
+ # Second case in default parameters and overwritted
+ if sec in overwrite and p in overwrite[sec]:
+ assert overwrite[sec][p] == params[p]
+ continue
+
+ # Third case in default
+ assert params[p] == default[sec][p]
+
+print("Checking if all the elements exist")
+# Finally ensure that all the elements are present
+for sec in default.keys():
+ params = default[sec]
+
+ for p in params.keys():
+ if sec not in over:
+ raise Exception("Section {} not found".format(sec))
+
+ if p not in over[sec]:
+ raise Exception("Parameter {}:{} not found".format(sec, p))
diff --git a/tests/run_tests.sh b/tests/run_tests.sh
new file mode 100644
index 0000000000000000000000000000000000000000..6d458ebc459d3846e444c5465270e172f074f06a
--- /dev/null
+++ b/tests/run_tests.sh
@@ -0,0 +1,9 @@
+#!/bin/bash
+
+set -e
+
+
+for filename in ./*.py;
+do
+ python $filename
+done
diff --git a/tests/test_cooling.py b/tests/test_cooling.py
new file mode 100644
index 0000000000000000000000000000000000000000..64c489befb53c8f6eb912ef194cdc865bdffc24b
--- /dev/null
+++ b/tests/test_cooling.py
@@ -0,0 +1,91 @@
+#!/usr/bin/env python3
+# ########################################################################
+# This file is part of PYSWIFT.
+# Copyright (c) 2018 Loic Hausammann (loic.hausammann@epfl.ch)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+# ########################################################################
+
+from pyswiftsim.libcooling import coolingRate
+import pyswiftsim
+import numpy as np
+from astropy import units, constants
+import yaml
+
+# define parameters
+filename = "test_cooling.yml"
+# adiabatic index
+gamma = 5. / 3.
+
+# download cooling table
+pyswiftsim.downloadCoolingTable()
+
+
+def computeCooling(filename):
+
+ # Read parameter file
+ f = open(filename, "r")
+ data = yaml.load(f)
+ f.close()
+
+ # Get the units
+ swift_units = data["InternalUnitSystem"]
+ # length unit
+ u_l = float(swift_units["UnitLength_in_cgs"])
+ # velocity unit
+ u_v = float(swift_units["UnitVelocity_in_cgs"])
+ # mass unit
+ u_m = float(swift_units["UnitMass_in_cgs"])
+ # time unit
+ u_t = u_l / u_v
+ # energy unit
+ u_e = u_m * u_l**2 / u_t**2
+
+ # boltzman constant
+ k_B = constants.k_B.to("erg / K") * units.K / (units.erg * u_e)
+ k_B = float(k_B)
+ # proton mass
+ m_p = constants.m_p.to("g") / (units.g * u_m)
+ m_p = float(m_p)
+
+ def internalEnergy(T):
+ u = (k_B / m_p) * T
+ u /= (gamma - 1.)
+ return u
+
+ # generates data
+ # time step in Myr
+ dt = 1e-3
+ dt = float(units.Myr.to('s')) / u_t
+
+ # density in atom / cm3
+ density = np.array([1.], dtype=np.float32)
+ density *= m_p * u_l**3
+
+ # Temperature in K
+ T = np.array([1e4], dtype=np.float32)
+ u = internalEnergy(T)
+
+ # compute the cooling rate
+ with_cosmo = False
+ rate = coolingRate(filename, density, u, with_cosmo, dt)
+
+ # go back to cgs
+ rate /= u_e / u_t
+
+ print("Cooling done: {} erg/s".format(rate))
+
+
+with pyswiftsim.ParameterFile() as filename:
+ computeCooling(filename)
diff --git a/tests/test_initial_mass_function.py b/tests/test_initial_mass_function.py
new file mode 100644
index 0000000000000000000000000000000000000000..ca80c505bb39844035953934cdb84bf221591068
--- /dev/null
+++ b/tests/test_initial_mass_function.py
@@ -0,0 +1,47 @@
+#!/usr/bin/env python3
+# ########################################################################
+# This file is part of PYSWIFT.
+# Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+# ########################################################################
+
+from pyswiftsim import libstellar
+import pyswiftsim
+
+import numpy as np
+import matplotlib.pyplot as plt
+plt.rc('text', usetex=True)
+
+
+N = 100
+
+
+if __name__ == "__main__":
+
+ with pyswiftsim.ParameterFile() as filename:
+
+ parser = pyswiftsim.parseYamlFile(filename)
+ imf = parser["GEARInitialMassFunction"]
+ mass_limits = imf["mass_limits"]
+ mass_min = float(mass_limits[0])
+ mass_max = float(mass_limits[1])
+
+ mass = np.linspace(mass_min, mass_max, N, dtype=np.float32)
+
+ print("Computing initial mass function...")
+
+ imf = libstellar.initialMassFunction(filename, mass)
+
+ print("IMF obtained: {}".format(imf))
diff --git a/tests/test_lifetime.py b/tests/test_lifetime.py
new file mode 100644
index 0000000000000000000000000000000000000000..9a1816e22438b5a857a9945e571961af280cccbc
--- /dev/null
+++ b/tests/test_lifetime.py
@@ -0,0 +1,57 @@
+#!/usr/bin/env python3
+# ########################################################################
+# This file is part of PYSWIFT.
+# Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+# ########################################################################
+
+from pyswiftsim import libstellar
+import pyswiftsim
+
+import numpy as np
+import matplotlib.pyplot as plt
+plt.rc('text', usetex=True)
+
+
+N = 100
+solar_metallicity = 0.02041
+
+
+if __name__ == "__main__":
+ with pyswiftsim.ParameterFile() as filename:
+
+ parser = pyswiftsim.parseYamlFile(filename)
+ # Get masses
+ imf = parser["GEARInitialMassFunction"]
+ mass_limits = imf["mass_limits"]
+ mass_min = float(mass_limits[0])
+ mass_max = float(mass_limits[1])
+
+ mass = np.linspace(mass_min, mass_max, N, dtype=np.float32)
+
+ # get metallicity
+ Z = solar_metallicity * np.ones(N, dtype=np.float32)
+
+ print("Computing lifetime...")
+
+ lifetime = libstellar.lifetimeFromMass(filename, mass, Z)
+ mass_from_lifetime = libstellar.massFromLifetime(filename, lifetime, Z)
+
+ print("Mass used: {}".format(mass))
+
+ print("Lifetime obtained: {}".format(lifetime))
+
+ print("Relative error on the mass from the lifetime: {}".format(
+ (mass_from_lifetime - mass) / mass))