Skip to content
Snippets Groups Projects
Commit c5d0d2ba authored by Loic Hausammann's avatar Loic Hausammann
Browse files

Update code with changes in swift

parent 58287259
No related branches found
No related tags found
1 merge request!10Update SWIFT calls
Hide whitespace changes
Inline Side-by-side
examples/cooling_box/plot_cooling.py
+ 1
− 1
View file @ c5d0d2ba
......@@ -131,7 +131,7 @@ if __name__ == "__main__":
overwrite = {
"GrackleCooling": {
"WithMetalCooling": with_metals,
"with_metal_cooling": with_metals,
}
}
with pyswiftsim.ParameterFile(overwrite) as filename:
......
examples/cooling_rate/plot_cooling.py
+ 111
− 52
View file @ c5d0d2ba
......@@ -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)
pyswiftsim/__init__.py
+ 16
− 15
View file @ c5d0d2ba
......@@ -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
},
}
"""
......
scripts/pyswift_lifetime
+ 3
− 3
View file @ c5d0d2ba
......@@ -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
......
src/cooling_wrapper.c
+ 11
− 13
View file @ c5d0d2ba
......@@ -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(&params, &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;
......
src/cooling_wrapper.h
+ 4
− 18
View file @ c5d0d2ba
......@@ -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(&params, &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(&params, &us, &pconst, &hydro_props, \
&cooling); \
#endif // __PYSWIFTSIM_COOLING_H__
......
src/stellar_evolution_wrapper.c
+ 13
− 10
View file @ c5d0d2ba
......@@ -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, &params, &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);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment