diff --git a/examples/cooling_box/plot_cooling.py b/examples/cooling_box/plot_cooling.py
index 18bbb26903a25e87667887dc5beeacc8f3feb976..0a3f2ee4c5262c4cdeb3dc22058ae7953af98bd1 100644
--- a/examples/cooling_box/plot_cooling.py
+++ b/examples/cooling_box/plot_cooling.py
@@ -131,7 +131,7 @@ if __name__ == "__main__":
overwrite = {
"GrackleCooling": {
- "WithMetalCooling": with_metals,
+ "with_metal_cooling": with_metals,
}
}
with pyswiftsim.ParameterFile(overwrite) as filename:
diff --git a/examples/cooling_rate/plot_cooling.py b/examples/cooling_rate/plot_cooling.py
index 884f0592da44acc0444e14e41ea63dcba67bf4eb..0fe4d8b42a5c7f99645b9957120f738a36817e32 100644
--- a/examples/cooling_rate/plot_cooling.py
+++ b/examples/cooling_rate/plot_cooling.py
@@ -20,6 +20,8 @@
from pyswiftsim.libcooling import coolingRate
import pyswiftsim
+import matplotlib
+
from copy import deepcopy
import numpy as np
import matplotlib.pyplot as plt
@@ -27,8 +29,20 @@ from astropy import units, constants
from time import strftime, gmtime
import sys
+SMALL_SIZE = 13
+MEDIUM_SIZE = 15
+BIGGER_SIZE = 18
+
+plt.rc('font', size=SMALL_SIZE) # controls default text sizes
+plt.rc('axes', titlesize=SMALL_SIZE) # fontsize of the axes title
+plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
+plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels
+plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels
+plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize
+plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
plt.rc('text', usetex=True)
+
pyswiftsim.downloadCoolingTable()
# PARAMETERS
@@ -40,6 +54,12 @@ with_cosmo = False
# include metals?
with_metals = 1
+metallicity = 0.02041
+
+# cooling table
+cooling_tables = [
+ "./CloudyData_UVB=HM2012.h5",
+]
# hydrogen mass fraction
h_mass_frac = 0.76
@@ -50,6 +70,7 @@ if N_rho == 1:
density = np.array([1.])
else:
density = np.logspace(-6, 4, N_rho)
+ plot_cooling_length = False
# temperature in K
N_T = 1000
@@ -113,7 +134,7 @@ def generate_initial_condition(us):
return d
-def plot_solution(rate, data, us):
+def plot_solution(rates, data, us):
print("Plotting solution")
energy = data["energy"]
rho = data["density"]
@@ -131,55 +152,82 @@ def plot_solution(rate, data, us):
energy *= u_v**2
- rate *= u_v**2 / u_t
+ rates = [rate * u_v**2 / u_t for rate in rates]
# lambda cooling
- lam = rate * rho
+ lambdas = [rate * rho for rate in rates]
# do plot
if N_rho == 1:
# plot Lambda vs T
plt.figure()
- plt.loglog(T, np.abs(lam),
- label="SWIFT")
+ for i, lam in enumerate(lambdas):
+ plt.loglog(T, np.abs(lam),
+ label="Table %i" % i)
plt.xlabel("Temperature [K]")
plt.ylabel("$\\Lambda$ [erg s$^{-1}$ cm$^{3}$]")
+ plt.legend()
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(gamma * (gamma-1.) * energy) * 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 )")
+ r = []
+ for rate in rates:
+ shape = [N_rho, N_T]
+ energy = np.reshape(energy, shape)
+ rate = np.reshape(rate, shape)
+
+ cooling_time = energy / rate
+ cooling_length = np.sqrt(gamma * (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)
+
+ levels = 500
+ if plot_cooling_length:
+ _min = -7
+ _max = 7
+ levels = np.linspace(_min, _max, levels)
+ z = cooling_length
+ else:
+ z = np.log10(np.abs(rate))
+
+ plt.figure()
+ r.append(rate)
+ plt.contourf(np.log10(rho), np.log10(T), z,
+ levels,
+ cmap=plt.get_cmap(colormap))
+ plt.xlabel("Log( Density [atom/cm3] )")
+ plt.ylabel("Log( Temperature [K] )")
+
+ cbar = plt.colorbar()
+ if plot_cooling_length:
+ tc = np.arange(_min, _max, 1.5)
+ cbar.set_ticks(tc)
+ cbar.set_ticklabels(tc)
+ cbar.ax.set_ylabel("Log( Cooling Length [kpc] )")
+ else:
+ cbar.ax.set_ylabel("Log( Rate [s / erg] )")
+
+ if len(r) == 2:
+ plt.figure()
+ plt.contourf(rho, T, np.log10(np.abs((r[1] - r[0]) / r[1])), 500,
+ 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()
date = strftime("%d %b %Y", gmtime())
plt.title("Date: {}".format(date))
@@ -188,24 +236,35 @@ def plot_solution(rate, data, us):
if __name__ == "__main__":
- overwrite = {
- "GrackleCooling": {
- "WithMetalCooling": with_metals,
+ rates = []
+
+ for cooling_table in cooling_tables:
+ overwrite = {
+ "GrackleCooling": {
+ "with_metal_cooling": with_metals,
+ "redshift": 0,
+ "cloudy_table": cooling_table,
+ },
+ "GEARChemistry": {
+ "initial_metallicity": metallicity,
+ },
}
- }
- with pyswiftsim.ParameterFile(overwrite) as filename:
- parser = pyswiftsim.parseYamlFile(filename)
- us = parser["InternalUnitSystem"]
+ 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"])
+ rates.append(rate)
- plot_solution(rate, d, us)
+ plot_solution(rates, d, us)
diff --git a/pyswiftsim/__init__.py b/pyswiftsim/__init__.py
index 32b6d181e0160b2239e508181044d860bc81f4b6..de585175385c46080183b649f8d07a63798abf10 100644
--- a/pyswiftsim/__init__.py
+++ b/pyswiftsim/__init__.py
@@ -159,26 +159,27 @@ class ParameterFile:
},
"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,
- "ThermalTime_Myr": 5,
+ "cloudy_table": "CloudyData_UVB=HM2012.h5",
+ "with_UV_background": 1,
+ "redshift": 0,
+ "with_metal_cooling": 1,
+ "provide_volumetric_heating_rates": 0,
+ "provide_specific_heating_rates": 0,
+ "self_shielding_method": 0,
+ "convergence_limit": 1e-2,
+ "max_steps": 20000,
+ "omega": 0.8,
+ "thermal_time_myr": 5,
},
- "GearChemistry": {
- "InitialMetallicity": 0.01295,
+ "GEARChemistry": {
+ "initial_metallicity": 0.01295,
},
"GEARFeedback": {
- "SupernovaeEnergy_erg": 1e51,
- "YieldsTable": "chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5",
+ "supernovae_energy_erg": 1e51,
+ "yields_table": "chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5",
+ "discrete_yields": 0
},
}
"""
diff --git a/scripts/pyswift_lifetime b/scripts/pyswift_lifetime
index b576e86d91b2dd1dc1643ac89af9dd368a2951ad..da0a88dfd0fae15762f91f070ab0639813e767d2 100644
--- a/scripts/pyswift_lifetime
+++ b/scripts/pyswift_lifetime
@@ -92,9 +92,9 @@ if __name__ == "__main__":
# transform into internal units
units = parser["InternalUnitSystem"]
- unit_mass = units["UnitMass_in_cgs"]
- unit_vel = units["UnitVelocity_in_cgs"]
- unit_length = units["UnitLength_in_cgs"]
+ unit_mass = float(units["UnitMass_in_cgs"])
+ unit_vel = float(units["UnitVelocity_in_cgs"])
+ unit_length = float(units["UnitLength_in_cgs"])
unit_time = unit_length / unit_vel
mass *= solMass_in_cgs / unit_mass
diff --git a/src/cooling_wrapper.c b/src/cooling_wrapper.c
index 116d94159e7a4fae312a72a3403a374b9bf394e5..dbe718c9f17cf5c7f5ef01f313b5d694bd60ca0c 100644
--- a/src/cooling_wrapper.c
+++ b/src/cooling_wrapper.c
@@ -19,11 +19,6 @@
#include "cooling_wrapper.h"
#include "pyswiftsim_tools.h"
-/* TODO Remove this */
-struct cooling_function_data cooling;
-int initialized;
-
-
/**
* @brief Set the cooling element fractions
*
@@ -79,7 +74,6 @@ void pycooling_set_fractions(struct xpart *xp, PyArrayObject* frac, const int id
*/
PyObject* pycooling_rate(PyObject* self, PyObject* args) {
import_array();
-
char *filename;
@@ -125,7 +119,6 @@ PyObject* pycooling_rate(PyObject* self, PyObject* args) {
struct chemistry_global_data chemistry;
chemistry_init_backend(¶ms, &us, &pconst, &chemistry);
-
/* Init entropy floor */
struct entropy_floor_properties floor_props;
entropy_floor_init(&floor_props, &pconst, &us,
@@ -135,7 +128,7 @@ PyObject* pycooling_rate(PyObject* self, PyObject* args) {
struct xpart xp;
hydro_first_init_part(&p, &xp);
-
+
/* return object */
PyArrayObject *rate = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(energy), PyArray_DIMS(energy), NPY_FLOAT);
@@ -150,16 +143,20 @@ PyObject* pycooling_rate(PyObject* self, PyObject* args) {
chemistry_first_init_part(&pconst, &us, &cosmo, &chemistry, &p, &xp);
if (fractions != NULL)
- pycooling_set_fractions(&xp, fractions, i);
+ pycooling_set_fractions(&xp, fractions, i);
else {
- /* Need to set the mass in order to estimate the element fractions */
- hydro_set_mass(&p, p.rho);
- p.h = 1;
+ /* Need to set the mass in order to estimate the element fractions */
+ hydro_set_mass(&p, p.rho);
+ p.h = 1;
- cooling_first_init_part(&pconst, &us, &cosmo, &cooling, &p, &xp);
+ cooling_first_init_part(&pconst, &us, &hydro_props, &cosmo, &cooling, &p, &xp);
}
hydro_set_physical_internal_energy_dt(&p, &cosmo, 0);
+ for(int j = 0; j < GEAR_CHEMISTRY_ELEMENT_COUNT; j++) {
+ p.chemistry_data.smoothed_metal_mass_fraction[j] = p.chemistry_data.metal_mass_fraction[j];
+ }
+
/* compute cooling rate */
cooling_cool_part(&pconst, &us, &cosmo, &hydro_props,
@@ -171,6 +168,7 @@ PyObject* pycooling_rate(PyObject* self, PyObject* args) {
/* Cleanup */
cosmology_clean(&cosmo);
+ cooling_clean(&cooling);
return (PyObject *) rate;
diff --git a/src/cooling_wrapper.h b/src/cooling_wrapper.h
index e0b97a9c9f70d9033dfaac54ae75bb061f8dfb35..8491ff58ae39a044d0fa2e0f979f56ccb3a96481 100644
--- a/src/cooling_wrapper.h
+++ b/src/cooling_wrapper.h
@@ -21,26 +21,12 @@
#include "pyswiftsim_includes.h"
-/* A few global variables */
-extern struct cooling_function_data cooling;
-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; \
- } \
- \
+#define PYSWIFTSIM_INIT_COOLING() \
+ struct cooling_function_data cooling; \
+ cooling_init_backend(¶ms, &us, &pconst, &hydro_props, \
+ &cooling); \
#endif // __PYSWIFTSIM_COOLING_H__
diff --git a/src/stellar_evolution_wrapper.c b/src/stellar_evolution_wrapper.c
index 007e369a8a5550e95d3d15cdfbb0ab2fb882d22d..3a23c8d0fd80f6fbff27cc02c2b06e8fc267cf28 100644
--- a/src/stellar_evolution_wrapper.c
+++ b/src/stellar_evolution_wrapper.c
@@ -72,6 +72,9 @@ PyObject* pyinitial_mass_function(PyObject* self, PyObject* args) {
}
*imf_i = initial_mass_function_get_imf(&fp.stellar_model.imf, m);
+
+ /* change from internal units to solar mass */
+ *imf_i /= pconst.const_solar_mass;
}
return (PyObject *) imf;
@@ -280,7 +283,7 @@ PyObject* pysupernovae_ia_rate(PyObject* self, PyObject* args) {
m2 = pow(10, m2);
/* Need to reverse 1 and 2 due to the lifetime being a decreasing function */
- *r = supernovae_ia_get_number(&fp.stellar_model.snia, m2, m1) / (t2 - t1);
+ *r = supernovae_ia_get_number_per_unit_mass(&fp.stellar_model.snia, m2, m1) / (t2 - t1);
/* change units 1 / (solMass Myr) -> internal units */
*r /= pconst.const_solar_mass * pconst.const_year * 1e6;
@@ -360,7 +363,7 @@ PyObject* pysupernovae_ii_rate(PyObject* self, PyObject* args) {
m2 = pow(10, m2);
/* Need to reverse 1 and 2 due to the lifetime being a decreasing function */
- *r = supernovae_ii_get_number(&fp.stellar_model.snii, m2, m1) / (t2 - t1);
+ *r = supernovae_ii_get_number_per_unit_mass(&fp.stellar_model.snii, m2, m1) / (t2 - t1);
/* change units 1 / (solMass Myr) -> internal units */
*r /= pconst.const_solar_mass * pconst.const_year * 1e6;
@@ -421,11 +424,11 @@ PyObject* pysupernovae_ii_yields(PyObject* self, PyObject* args) {
/* return object */
const int ndims = 2;
- npy_intp dims[2] = {N, CHEMISTRY_ELEMENT_COUNT};
+ npy_intp dims[2] = {N, GEAR_CHEMISTRY_ELEMENT_COUNT};
PyArrayObject *yields_out = (PyArrayObject *)PyArray_SimpleNew(ndims, dims, NPY_FLOAT);
for(size_t i = 0; i < N; i++) {
- float yields[CHEMISTRY_ELEMENT_COUNT] = {0.};
+ float yields[GEAR_CHEMISTRY_ELEMENT_COUNT] = {0.};
const float m1 = *(float*) PyArray_GETPTR1(mass1, i) / pconst.const_solar_mass;
const float log_m1 = log10(m1);
const float m2 = *(float*) PyArray_GETPTR1(mass2, i) / pconst.const_solar_mass;
@@ -434,7 +437,7 @@ PyObject* pysupernovae_ii_yields(PyObject* self, PyObject* args) {
/* Compute the yields */
supernovae_ii_get_yields(&fp.stellar_model.snii, log_m1, log_m2, yields);
- for(int j = 0; j < CHEMISTRY_ELEMENT_COUNT; j++) {
+ for(int j = 0; j < GEAR_CHEMISTRY_ELEMENT_COUNT; j++) {
float *y = (float*) PyArray_GETPTR2(yields_out, i, j);
*y = yields[j];
}
@@ -608,12 +611,12 @@ PyObject* pysupernovae_ia_yields(PyObject* self, PyObject* args) {
/* return object */
const int ndims = 1;
- npy_intp dims[1] = {CHEMISTRY_ELEMENT_COUNT};
+ npy_intp dims[1] = {GEAR_CHEMISTRY_ELEMENT_COUNT};
PyArrayObject *yields_out = (PyArrayObject *)PyArray_SimpleNew(ndims, dims, NPY_FLOAT);
const float *yields = supernovae_ia_get_yields(&fp.stellar_model.snia);
- for(size_t i = 0; i < CHEMISTRY_ELEMENT_COUNT; i++) {
+ for(size_t i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
float *y = (float*) PyArray_GETPTR1(yields_out, i);
*y = yields[i] / fp.stellar_model.snia.mass_white_dwarf;
@@ -658,7 +661,7 @@ PyObject* pysupernovae_ia_mass_ejected(PyObject* self, PyObject* args) {
float *y = (float*) PyArray_GETPTR1(mass_ejected, 0);
- *y = supernovae_ia_get_ejected_mass_fraction_processed(&fp.stellar_model.snia) / fp.stellar_model.snia.mass_white_dwarf;
+ *y = supernovae_ia_get_ejected_mass_processed(&fp.stellar_model.snia) / fp.stellar_model.snia.mass_white_dwarf;
return (PyObject *) mass_ejected;
@@ -693,9 +696,9 @@ PyObject* pystellar_get_elements(PyObject* self, PyObject* args) {
feedback_props_init(&fp, &pconst, &us, ¶ms, &hydro_props, &cosmo);
/* Save the variables */
- PyObject *t = PyList_New(CHEMISTRY_ELEMENT_COUNT);
+ PyObject *t = PyList_New(GEAR_CHEMISTRY_ELEMENT_COUNT);
- for(int i = 0; i < CHEMISTRY_ELEMENT_COUNT; i++) {
+ for(int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
const char *name = stellar_evolution_get_element_name(&fp.stellar_model, i);
PyObject *string = PyUnicode_FromString(name);
int test = PyList_SetItem(t, i, string);