diff --git a/examples/Cooling/ConstantCosmoTempEvolution/plot_thermal_history.py b/examples/Cooling/ConstantCosmoTempEvolution/plot_thermal_history.py
index 1494102531104b252e3edaa467920db7383ac6e6..1f9ae6557d1a55bd75d3440679ae42d653da3399 100644
--- a/examples/Cooling/ConstantCosmoTempEvolution/plot_thermal_history.py
+++ b/examples/Cooling/ConstantCosmoTempEvolution/plot_thermal_history.py
@@ -68,21 +68,21 @@ for snap in snap_list:
z = np.append(z, data.metadata.z)
# Convert gas temperatures and densities to right units
- data.gas.temperature.convert_to_cgs()
+ data.gas.temperatures.convert_to_cgs()
# Get mean and standard deviation of temperature
- T_mean.append(np.mean(data.gas.temperature) * data.gas.temperature.units)
- T_std.append(np.std(data.gas.temperature) * data.gas.temperature.units)
+ T_mean.append(np.mean(data.gas.temperatures) * data.gas.temperatures.units)
+ T_std.append(np.std(data.gas.temperatures) * data.gas.temperatures.units)
# Get mean and standard deviation of density
- rho_mean.append(np.mean(data.gas.density) * data.gas.density.units)
- rho_std.append(np.std(data.gas.density) * data.gas.density.units)
+ rho_mean.append(np.mean(data.gas.densities) * data.gas.densities.units)
+ rho_std.append(np.std(data.gas.densities) * data.gas.densities.units)
## Turn into numpy arrays
-T_mean = np.array(T_mean) * data.gas.temperature.units
-T_std = np.array(T_std) * data.gas.temperature.units
-rho_mean = np.array(rho_mean) * data.gas.density.units
-rho_std = np.array(rho_std) * data.gas.density.units
+T_mean = np.array(T_mean) * data.gas.temperatures.units
+T_std = np.array(T_std) * data.gas.temperatures.units
+rho_mean = np.array(rho_mean) * data.gas.densities.units
+rho_std = np.array(rho_std) * data.gas.densities.units
## Put Density into units of mean baryon density today
diff --git a/examples/Cooling/CoolingBox/plotEnergy.py b/examples/Cooling/CoolingBox/plotEnergy.py
index 9c7af57d3d9dfdcfa222e9f77701f230d25f9ddc..6234c319a3001fa4e364639bb2c5480a31cf99dd 100644
--- a/examples/Cooling/CoolingBox/plotEnergy.py
+++ b/examples/Cooling/CoolingBox/plotEnergy.py
@@ -87,7 +87,7 @@ temp_snap = np.zeros(N)
time_snap_cgs = np.zeros(N)
for i in range(N):
snap = File(files[i], 'r')
- u = snap["/PartType0/InternalEnergy"][:] * snap["/PartType0/Masses"][:]
+ u = snap["/PartType0/InternalEnergies"][:] * snap["/PartType0/Masses"][:]
u = sum(u) / total_mass[0]
temp_snap[i] = energyUnits(u)
time_snap_cgs[i] = snap["/Header"].attrs["Time"] * unit_time
diff --git a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py
index 9cde36cb05b88e60b3bf3527f3857a3774bf5dca..cb624be3cfb1f0f803cfce7a14fb1a772cccf515 100644
--- a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py
+++ b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py
@@ -94,15 +94,17 @@ def get_data(handle: float, n_snaps: int):
t0 = data.metadata.t.to(Myr).value
times.append(data.metadata.t.to(Myr).value - t0)
- temps.append(np.mean(data.gas.temperature.to(K).value))
+ temps.append(np.mean(data.gas.temperatures.to(K).value))
densities.append(
- np.mean(data.gas.density.to(mh / cm ** 3).value)
+ np.mean(data.gas.densities.to(mh / cm ** 3).value)
/ (data.metadata.scale_factor ** 3)
)
try:
energies.append(
- np.mean((data.gas.internal_energy * data.gas.masses).to(erg).value)
+ np.mean(
+ (data.gas.internal_energies * data.gas.masses).to(erg).value
+ )
* data.metadata.scale_factor ** (2)
)
radiated_energies.append(
diff --git a/examples/Cooling/FeedbackEvent_3D/plotSolution.py b/examples/Cooling/FeedbackEvent_3D/plotSolution.py
index fe6f93996dafee2b6e81a70d4786374d86355f6f..6631fff2244e73244d0311f4e823c42dfcea7609 100644
--- a/examples/Cooling/FeedbackEvent_3D/plotSolution.py
+++ b/examples/Cooling/FeedbackEvent_3D/plotSolution.py
@@ -50,7 +50,6 @@ except:
plt.rcParams.update(rcParams)
-
def get_data_dump(metadata):
"""
Gets a big data dump from the SWIFT metadata
@@ -104,10 +103,10 @@ r = r.to(kpc).value
data = dict(
x=r,
v_r=v_r,
- u=sim.gas.temperature.to(K).value,
- S=sim.gas.entropy.to(keV / K).value,
- P=sim.gas.pressure.to(kPa).value,
- rho=sim.gas.density.to(mh / (cm ** 3)).value,
+ u=sim.gas.temperatures.to(K).value,
+ S=sim.gas.entropies.to(keV / K).value,
+ P=sim.gas.pressures.to(kPa).value,
+ rho=sim.gas.densities.to(mh / (cm ** 3)).value,
)
# Try to add on the viscosity and diffusion.
@@ -156,8 +155,13 @@ log = dict(
v_r=False, v_phi=False, u=False, S=False, P=False, rho=False, visc=False, diff=False
)
ylim = dict(
- v_r=[-4, 25], u=[4750, 6000], rho=[0.09, 0.16], visc=[0, 2.0], diff=[0, 0.25],
- P=[3e-18, 10e-18], S=[-0.5e60, 4e60]
+ v_r=[-4, 25],
+ u=[4750, 6000],
+ rho=[0.09, 0.16],
+ visc=[0, 2.0],
+ diff=[0, 0.25],
+ P=[3e-18, 10e-18],
+ S=[-0.5e60, 4e60],
)
current_axis = 0
diff --git a/examples/Cosmology/ComovingSodShock_1D/plotSolution.py b/examples/Cosmology/ComovingSodShock_1D/plotSolution.py
index 95674c04bfafd0cd549b69814df82f9a4f80a949..b09873339fc4ded024c70f66301c5cce762b97fc 100644
--- a/examples/Cosmology/ComovingSodShock_1D/plotSolution.py
+++ b/examples/Cosmology/ComovingSodShock_1D/plotSolution.py
@@ -82,12 +82,12 @@ git = str(sim["Code"].attrs["Git Revision"])
x = sim["/PartType0/Coordinates"][:,0]
v = sim["/PartType0/Velocities"][:,0] * anow
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
try:
- alpha = sim["/PartType0/Viscosity"][:]
+ alpha = sim["/PartType0/ViscosityParameters"][:]
plot_alpha = True
except:
plot_alpha = False
diff --git a/examples/Cosmology/ComovingSodShock_2D/plotSolution.py b/examples/Cosmology/ComovingSodShock_2D/plotSolution.py
index 8adb3cf5c550ab9724f6a8f34c1a1260a25712e1..ec28b9449bfa6b00d3088952a6b1bf871462f86c 100644
--- a/examples/Cosmology/ComovingSodShock_2D/plotSolution.py
+++ b/examples/Cosmology/ComovingSodShock_2D/plotSolution.py
@@ -83,10 +83,10 @@ git = sim["Code"].attrs["Git Revision"]
x = sim["/PartType0/Coordinates"][:,0]
v = sim["/PartType0/Velocities"][:,0] * anow
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
N = 1000 # Number of points
x_min = -1.
diff --git a/examples/Cosmology/ComovingSodShock_3D/plotSolution.py b/examples/Cosmology/ComovingSodShock_3D/plotSolution.py
index d85f4be9a49d108d133928a81ea4482fa9099792..34abae364f5421a0436352b9564537f2f31e8742 100644
--- a/examples/Cosmology/ComovingSodShock_3D/plotSolution.py
+++ b/examples/Cosmology/ComovingSodShock_3D/plotSolution.py
@@ -83,19 +83,19 @@ git = sim["Code"].attrs["Git Revision"]
x = sim["/PartType0/Coordinates"][:,0]
v = sim["/PartType0/Velocities"][:,0] * anow
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
try:
- diffusion = sim["/PartType0/Diffusion"][:]
+ diffusion = sim["/PartType0/DiffusionParameters"][:]
plot_diffusion = True
except:
plot_diffusion = False
try:
- viscosity = sim["/PartType0/Viscosity"][:]
+ viscosity = sim["/PartType0/ViscosityParameters"][:]
plot_viscosity = True
except:
plot_viscosity = False
diff --git a/examples/Cosmology/ConstantCosmoVolume/plotSolution.py b/examples/Cosmology/ConstantCosmoVolume/plotSolution.py
index f77889d7cb19c230accf25290b88a321e0f59616..6df0bfa4e1f7c4932f652d82a8ca95f5c54b0e9f 100644
--- a/examples/Cosmology/ConstantCosmoVolume/plotSolution.py
+++ b/examples/Cosmology/ConstantCosmoVolume/plotSolution.py
@@ -109,19 +109,19 @@ for i in range(119):
z[i] = sim["/Cosmology"].attrs["Redshift"][0]
a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
- S = sim["/PartType0/Entropy"][:]
+ S = sim["/PartType0/Entropies"][:]
S_mean[i] = np.mean(S)
S_std[i] = np.std(S)
- u = sim["/PartType0/InternalEnergy"][:]
+ u = sim["/PartType0/InternalEnergies"][:]
u_mean[i] = np.mean(u)
u_std[i] = np.std(u)
- P = sim["/PartType0/Pressure"][:]
+ P = sim["/PartType0/Pressures"][:]
P_mean[i] = np.mean(P)
P_std[i] = np.std(P)
- rho = sim["/PartType0/Density"][:]
+ rho = sim["/PartType0/Densities"][:]
rho_mean[i] = np.mean(rho)
rho_std[i] = np.std(rho)
diff --git a/examples/Cosmology/ZeldovichPancake_3D/plotSolution.py b/examples/Cosmology/ZeldovichPancake_3D/plotSolution.py
index eef247fb761e75f8dde8e8abe84075efbd7cb46a..285ed3ae218549468c309e9109ef53f10a2f66ab 100644
--- a/examples/Cosmology/ZeldovichPancake_3D/plotSolution.py
+++ b/examples/Cosmology/ZeldovichPancake_3D/plotSolution.py
@@ -78,10 +78,10 @@ gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
x = sim["/PartType0/Coordinates"][:,0]
v = sim["/PartType0/Velocities"][:,0]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
m = sim["/PartType0/Masses"][:]
try:
phi = sim["/PartType0/Potential"][:]
@@ -98,8 +98,8 @@ if os.path.exists(filename_g):
sim_g = h5py.File(filename_g, "r")
x_g = sim_g["/PartType0/Coordinates"][:,0]
v_g = sim_g["/PartType0/Velocities"][:,0]
- u_g = sim_g["/PartType0/InternalEnergy"][:]
- rho_g = sim_g["/PartType0/Density"][:]
+ u_g = sim_g["/PartType0/InternalEnergies"][:]
+ rho_g = sim_g["/PartType0/Densities"][:]
phi_g = sim_g["/PartType0/Potential"][:]
a_g = sim_g["/Header"].attrs["Time"]
print("Gadget Scale-factor:", a_g, "redshift:", 1/a_g - 1.)
diff --git a/examples/HydroTests/BlobTest_3D/makeMovie.py b/examples/HydroTests/BlobTest_3D/makeMovie.py
index 9ae4a538e000fa7006f093b67d325ff433e97089..55acdaab01e22ecf8bae0dd24967a348ceabf34b 100644
--- a/examples/HydroTests/BlobTest_3D/makeMovie.py
+++ b/examples/HydroTests/BlobTest_3D/makeMovie.py
@@ -24,7 +24,7 @@ resolution = 1024
snapshot_name = "blob"
cmap = "Spectral_r"
text_args = dict(color="black")
-# plot = "pressure"
+# plot = "pressures"
# name = "Pressure $P$"
plot = "density"
name = "Fluid Density $\\rho$"
diff --git a/examples/HydroTests/Diffusion_1D/plotSolution.py b/examples/HydroTests/Diffusion_1D/plotSolution.py
index 66c8ffc6418f06589a2918ae4d8ed460b0081972..a394e919b4e2e80728c97499c3c3544256a6ad2c 100644
--- a/examples/HydroTests/Diffusion_1D/plotSolution.py
+++ b/examples/HydroTests/Diffusion_1D/plotSolution.py
@@ -24,6 +24,7 @@ import numpy as np
try:
from scipy.integrate import solve_ivp
+
solve_ode = True
except:
solve_ode = False
@@ -32,6 +33,7 @@ from swiftsimio import load
matplotlib.use("Agg")
+
def solve_analytic(u_0, u_1, t_0, t_1, alpha=0.1):
"""
Solves the analytic equation:
@@ -63,7 +65,12 @@ def solve_analytic(u_0, u_1, t_0, t_1, alpha=0.1):
return np.array([-1.0 * common, 1.0 * common])
- ret = solve_ivp(gradient, t_span=[t_0.value, t_1.value], y0=[u_0.value, u_1.value], t_eval=np.linspace(t_0.value, t_1.value, 100))
+ ret = solve_ivp(
+ gradient,
+ t_span=[t_0.value, t_1.value],
+ y0=[u_0.value, u_1.value],
+ t_eval=np.linspace(t_0.value, t_1.value, 100),
+ )
t = ret.t
high = ret.y[1]
@@ -81,7 +88,7 @@ def get_data_dump(metadata):
viscosity = metadata.viscosity_info
except:
viscosity = "No info"
-
+
try:
diffusion = metadata.diffusion_info
except:
@@ -136,6 +143,7 @@ def setup_axes(size=[8, 8], dpi=300):
return fig, ax
+
def mean_std_max_min(data):
"""
Returns:
@@ -157,7 +165,7 @@ def extract_plottables_u(data_list):
"""
data = [
- np.diff(x.gas.internal_energy.value) / np.mean(x.gas.internal_energy.value)
+ np.diff(x.gas.internal_energies.value) / np.mean(x.gas.internal_energies.value)
for x in data_list
]
@@ -175,8 +183,10 @@ def extract_plottables_x(data_list):
dx = boxsize / n_part
original_x = np.arange(n_part, dtype=float) * dx + (0.5 * dx)
-
- deviations = [1e6 * abs(original_x - x.gas.coordinates.value[:, 0]) / dx for x in data_list]
+
+ deviations = [
+ 1e6 * abs(original_x - x.gas.coordinates.value[:, 0]) / dx for x in data_list
+ ]
return mean_std_max_min(deviations)
@@ -187,7 +197,7 @@ def extract_plottables_rho(data_list):
mean, stdev, max, min * 1e6 deviations from mean density
"""
- P = [x.gas.density.value for x in data_list]
+ P = [x.gas.densities.value for x in data_list]
mean_P = [np.mean(x) for x in P]
deviations = [1e6 * (x - y) / x for x, y in zip(mean_P, P)]
@@ -241,28 +251,34 @@ def make_plot(start: int, stop: int, handle: str):
if solve_ode:
times_ode, diff = solve_analytic(
- u_0=data_list[0].gas.internal_energy.min(),
- u_1=data_list[0].gas.internal_energy.max(),
+ u_0=data_list[0].gas.internal_energies.min(),
+ u_1=data_list[0].gas.internal_energies.max(),
t_0=t[0],
t_1=t[-1],
alpha=(
- np.sqrt(5.0/3.0 * (5.0/3.0 - 1.0)) *
- alpha / data_list[0].gas.smoothing_length[0].value
- )
+ np.sqrt(5.0 / 3.0 * (5.0 / 3.0 - 1.0))
+ * alpha
+ / data_list[0].gas.smoothing_length[0].value
+ ),
)
ax[1].plot(
times_ode,
- (diff) / np.mean(data_list[0].gas.internal_energy),
+ (diff) / np.mean(data_list[0].gas.internal_energies),
label="Analytic",
linestyle="dotted",
- c="C3"
+ c="C3",
)
- #import pdb;pdb.set_trace()
+ # import pdb;pdb.set_trace()
ax[2].fill_between(
- t, x_means - x_stdevs, x_means + x_stdevs, color="C0", alpha=0.5, edgecolor="none"
+ t,
+ x_means - x_stdevs,
+ x_means + x_stdevs,
+ color="C0",
+ alpha=0.5,
+ edgecolor="none",
)
ax[2].plot(t, x_means, label="Mean", c="C0")
ax[2].plot(t, x_maxs, label="Max", linestyle="dashed", c="C1")
@@ -270,11 +286,18 @@ def make_plot(start: int, stop: int, handle: str):
try:
# Give diffusion info a go; this may not be present
- diff_means, diff_stdevs, diff_maxs, diff_mins = extract_plottables_diff(data_list)
+ diff_means, diff_stdevs, diff_maxs, diff_mins = extract_plottables_diff(
+ data_list
+ )
ax[3].set_ylabel(r"Diffusion parameter $\alpha_{diff}$")
ax[3].fill_between(
- t, diff_means - diff_stdevs, diff_means + diff_stdevs, color="C0", alpha=0.5, edgecolor="none"
+ t,
+ diff_means - diff_stdevs,
+ diff_means + diff_stdevs,
+ color="C0",
+ alpha=0.5,
+ edgecolor="none",
)
ax[3].plot(t, diff_means, label="Mean", c="C0")
ax[3].plot(t, diff_maxs, label="Max", linestyle="dashed", c="C1")
@@ -284,9 +307,16 @@ def make_plot(start: int, stop: int, handle: str):
# Diffusion info must not be present.
rho_means, rho_stdevs, rho_maxs, rho_mins = extract_plottables_rho(data_list)
- ax[3].set_ylabel("Deviation from mean density $(\\rho_i - \\bar{\\rho}) / \\bar{\\rho}$ [$\\times 10^{6}$]")
+ ax[3].set_ylabel(
+ "Deviation from mean density $(\\rho_i - \\bar{\\rho}) / \\bar{\\rho}$ [$\\times 10^{6}$]"
+ )
ax[3].fill_between(
- t, rho_means - rho_stdevs, rho_means + rho_stdevs, color="C0", alpha=0.5, edgecolor="none"
+ t,
+ rho_means - rho_stdevs,
+ rho_means + rho_stdevs,
+ color="C0",
+ alpha=0.5,
+ edgecolor="none",
)
ax[3].plot(t, rho_means, label="Mean", c="C0")
ax[3].plot(t, rho_maxs, label="Max", linestyle="dashed", c="C1")
diff --git a/examples/HydroTests/EvrardCollapse_3D/plotSolution.py b/examples/HydroTests/EvrardCollapse_3D/plotSolution.py
index 8422b9c45fd573f3d0ae36324d6e39ab23cceb25..b405771a4792e5ca4fef181b3787952bf0078d67 100644
--- a/examples/HydroTests/EvrardCollapse_3D/plotSolution.py
+++ b/examples/HydroTests/EvrardCollapse_3D/plotSolution.py
@@ -75,10 +75,10 @@ x = sqrt((coords[:,0] - 0.5 * boxSize)**2 + (coords[:,1] - 0.5 * boxSize)**2 + \
(coords[:,2] - 0.5 * boxSize)**2)
vels = sim["/PartType0/Velocities"]
v = sqrt(vels[:,0]**2 + vels[:,1]**2 + vels[:,2]**2)
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
# Bin the data
x_bin_edge = logspace(-3., log10(2.), 100)
diff --git a/examples/HydroTests/Gradients/plot.py b/examples/HydroTests/Gradients/plot.py
index d6750ffc581437ebbf402ec44bcb1d14cb82a698..7b4248c9fc9c9d6b9c5b15410185d6489849bed8 100644
--- a/examples/HydroTests/Gradients/plot.py
+++ b/examples/HydroTests/Gradients/plot.py
@@ -30,7 +30,7 @@ inputfile = sys.argv[1]
outputfile = "gradients_{0}.png".format(sys.argv[2])
f = h5py.File(inputfile, "r")
-rho = np.array(f["/PartType0/Density"])
+rho = np.array(f["/PartType0/Densities"])
gradrho = np.array(f["/PartType0/GradDensity"])
coords = np.array(f["/PartType0/Coordinates"])
diff --git a/examples/HydroTests/GreshoVortex_2D/plotSolution.py b/examples/HydroTests/GreshoVortex_2D/plotSolution.py
index d497a6b297bf38b39cf85a9107a769c20f815b77..2d4697b6ffaac0639da67ee90d824c75791ea573 100644
--- a/examples/HydroTests/GreshoVortex_2D/plotSolution.py
+++ b/examples/HydroTests/GreshoVortex_2D/plotSolution.py
@@ -100,10 +100,10 @@ r = sqrt(x**2 + y**2)
v_r = (x * vel[:,0] + y * vel[:,1]) / r
v_phi = (-y * vel[:,0] + x * vel[:,1]) / r
v_norm = sqrt(vel[:,0]**2 + vel[:,1]**2)
-rho = sim["/PartType0/Density"][:]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Densities"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
# Bin te data
r_bin_edge = np.arange(0., 1., 0.02)
diff --git a/examples/HydroTests/GreshoVortex_3D/plotSolution.py b/examples/HydroTests/GreshoVortex_3D/plotSolution.py
index 545440c997d9ebc3ab11562d0a7d9fa143e23ed2..20beab7514759c764f5ca7c379183506b764a819 100644
--- a/examples/HydroTests/GreshoVortex_3D/plotSolution.py
+++ b/examples/HydroTests/GreshoVortex_3D/plotSolution.py
@@ -103,19 +103,19 @@ r = sqrt(x**2 + y**2)
v_r = (x * vel[:,0] + y * vel[:,1]) / r
v_phi = (-y * vel[:,0] + x * vel[:,1]) / r
v_norm = sqrt(vel[:,0]**2 + vel[:,1]**2)
-rho = sim["/PartType0/Density"][:]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Densities"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
try:
- diffusion = sim["/PartType0/Diffusion"][:]
+ diffusion = sim["/PartType0/DiffusionParameters"][:]
plot_diffusion = True
except:
plot_diffusion = False
try:
- viscosity = sim["/PartType0/Viscosity"][:]
+ viscosity = sim["/PartType0/ViscosityParameters"][:]
plot_viscosity = True
except:
plot_viscosity = False
diff --git a/examples/HydroTests/InteractingBlastWaves_1D/plotSolution.py b/examples/HydroTests/InteractingBlastWaves_1D/plotSolution.py
index 1719162dec34626d6f4ecb8267c4d06f85b3db26..d617fb239ce21acab73b5cb057dd3cdf4b260d59 100644
--- a/examples/HydroTests/InteractingBlastWaves_1D/plotSolution.py
+++ b/examples/HydroTests/InteractingBlastWaves_1D/plotSolution.py
@@ -55,11 +55,11 @@ snap = int(sys.argv[1])
# Open the file and read the relevant data
file = h5py.File("interactingBlastWaves_{0:04d}.hdf5".format(snap), "r")
x = file["/PartType0/Coordinates"][:,0]
-rho = file["/PartType0/Density"]
+rho = file["/PartType0/Densities"]
v = file["/PartType0/Velocities"][:,0]
-u = file["/PartType0/InternalEnergy"]
-S = file["/PartType0/Entropy"]
-P = file["/PartType0/Pressure"]
+u = file["/PartType0/InternalEnergies"]
+S = file["/PartType0/Entropies"]
+P = file["/PartType0/Pressures"]
time = file["/Header"].attrs["Time"][0]
scheme = file["/HydroScheme"].attrs["Scheme"]
diff --git a/examples/HydroTests/KelvinHelmholtz_2D/plotSolution.py b/examples/HydroTests/KelvinHelmholtz_2D/plotSolution.py
index 77ab6fb244da25d13760f90653fac7eac11a0ee7..f599fcb784633b2d6765ea79767fc658196faa5f 100644
--- a/examples/HydroTests/KelvinHelmholtz_2D/plotSolution.py
+++ b/examples/HydroTests/KelvinHelmholtz_2D/plotSolution.py
@@ -77,10 +77,10 @@ x = pos[:,0] - boxSize / 2
y = pos[:,1] - boxSize / 2
vel = sim["/PartType0/Velocities"][:,:]
v_norm = sqrt(vel[:,0]**2 + vel[:,1]**2)
-rho = sim["/PartType0/Density"][:]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Densities"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
# Plot the interesting quantities
figure()
diff --git a/examples/HydroTests/Noh_1D/plotSolution.py b/examples/HydroTests/Noh_1D/plotSolution.py
index 25b9b2f16b24cba5def592a5cf00dbae82195ef7..7f0b5d403ef816b0dda57823010472476a7ecc32 100644
--- a/examples/HydroTests/Noh_1D/plotSolution.py
+++ b/examples/HydroTests/Noh_1D/plotSolution.py
@@ -69,10 +69,10 @@ git = sim["Code"].attrs["Git Revision"]
x = sim["/PartType0/Coordinates"][:,0]
v = sim["/PartType0/Velocities"][:,0]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
N = 1001 # Number of points
x_min = -1
diff --git a/examples/HydroTests/Noh_2D/plotSolution.py b/examples/HydroTests/Noh_2D/plotSolution.py
index 775ddf4e8a7954c14034ad51a6b66622c41a6996..b53212c4688ec790ad8f3f83f81243f9ec52266d 100644
--- a/examples/HydroTests/Noh_2D/plotSolution.py
+++ b/examples/HydroTests/Noh_2D/plotSolution.py
@@ -71,10 +71,10 @@ x = sim["/PartType0/Coordinates"][:,0]
y = sim["/PartType0/Coordinates"][:,1]
vx = sim["/PartType0/Velocities"][:,0]
vy = sim["/PartType0/Velocities"][:,1]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
r = np.sqrt((x-1)**2 + (y-1)**2)
v = -np.sqrt(vx**2 + vy**2)
diff --git a/examples/HydroTests/Noh_3D/plotSolution.py b/examples/HydroTests/Noh_3D/plotSolution.py
index 386b9f728b5e8d8e38fb7ec9aeaa336d194e35dd..20e8ca805de1cb700b8b462ae27495080f5d3268 100644
--- a/examples/HydroTests/Noh_3D/plotSolution.py
+++ b/examples/HydroTests/Noh_3D/plotSolution.py
@@ -74,10 +74,10 @@ z = sim["/PartType0/Coordinates"][:,2]
vx = sim["/PartType0/Velocities"][:,0]
vy = sim["/PartType0/Velocities"][:,1]
vz = sim["/PartType0/Velocities"][:,2]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
r = np.sqrt((x-1)**2 + (y-1)**2 + (z-1)**2)
v = -np.sqrt(vx**2 + vy**2 + vz**2)
diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/plotInitialProfile.py b/examples/HydroTests/Rayleigh-Taylor_2D/plotInitialProfile.py
index e89c4e8525fe4c88e517acbd453b0941f8f573c8..8928a6e597adbe4e52f905ba95fa68f424b6cabb 100644
--- a/examples/HydroTests/Rayleigh-Taylor_2D/plotInitialProfile.py
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/plotInitialProfile.py
@@ -12,8 +12,8 @@ f = load(filename)
# Get data from snapshot
x, y, _ = f.gas.coordinates.value.T
-rho = f.gas.density.value
-a = f.gas.entropy.value
+rho = f.gas.densities.value
+a = f.gas.entropies.value
# Get analytical solution
y_an = np.linspace(0, makeIC.box_size[1], N)
diff --git a/examples/HydroTests/SedovBlast_1D/plotSolution.py b/examples/HydroTests/SedovBlast_1D/plotSolution.py
index c6d4a989da7493f7b500946610eea8832696bf6f..d82ad9e94610a916d900ea93863b2881757e73b3 100644
--- a/examples/HydroTests/SedovBlast_1D/plotSolution.py
+++ b/examples/HydroTests/SedovBlast_1D/plotSolution.py
@@ -78,13 +78,13 @@ x = pos[:,0] - boxSize / 2
vel = sim["/PartType0/Velocities"][:,:]
r = abs(x)
v_r = x * vel[:,0] / r
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
try:
- alpha = sim["/PartType0/Viscosity"][:]
+ alpha = sim["/PartType0/ViscosityParameters"][:]
plot_alpha = True
except:
plot_alpha = False
diff --git a/examples/HydroTests/SedovBlast_2D/plotSolution.py b/examples/HydroTests/SedovBlast_2D/plotSolution.py
index 2b5de6f32b8673bbc825fbb5236f4e2ab3b4f408..6f504a09c9432368ce141ec0d28c28699f5ba7f3 100644
--- a/examples/HydroTests/SedovBlast_2D/plotSolution.py
+++ b/examples/HydroTests/SedovBlast_2D/plotSolution.py
@@ -80,10 +80,10 @@ y = pos[:,1] - boxSize / 2
vel = sim["/PartType0/Velocities"][:,:]
r = sqrt(x**2 + y**2)
v_r = (x * vel[:,0] + y * vel[:,1]) / r
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
# Bin te data
r_bin_edge = np.arange(0., 0.5, 0.01)
diff --git a/examples/HydroTests/SedovBlast_3D/plotSolution.py b/examples/HydroTests/SedovBlast_3D/plotSolution.py
index b0f2e08441b3fa550e61602ba852228a04362fbc..fec4f1101406be5803a3f1601812d1cb85275409 100644
--- a/examples/HydroTests/SedovBlast_3D/plotSolution.py
+++ b/examples/HydroTests/SedovBlast_3D/plotSolution.py
@@ -81,19 +81,19 @@ z = pos[:,2] - boxSize / 2
vel = sim["/PartType0/Velocities"][:,:]
r = sqrt(x**2 + y**2 + z**2)
v_r = (x * vel[:,0] + y * vel[:,1] + z * vel[:,2]) / r
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
try:
- diffusion = sim["/PartType0/Diffusion"][:]
+ diffusion = sim["/PartType0/DiffusionParameters"][:]
plot_diffusion = True
except:
plot_diffusion = False
try:
- viscosity = sim["/PartType0/Viscosity"][:]
+ viscosity = sim["/PartType0/ViscosityParameters"][:]
plot_viscosity = True
except:
plot_viscosity = False
diff --git a/examples/HydroTests/SineWavePotential_1D/plotSolution.py b/examples/HydroTests/SineWavePotential_1D/plotSolution.py
index 3bb889aaabd3cdac0274afb09647d0e3aebb02cc..ae99d98aaaa11d9c68473b74106054963a075895 100644
--- a/examples/HydroTests/SineWavePotential_1D/plotSolution.py
+++ b/examples/HydroTests/SineWavePotential_1D/plotSolution.py
@@ -43,9 +43,9 @@ fileName = sys.argv[1]
file = h5py.File(fileName, 'r')
coords = np.array(file["/PartType0/Coordinates"])
-rho = np.array(file["/PartType0/Density"])
-P = np.array(file["/PartType0/Pressure"])
-u = np.array(file["/PartType0/InternalEnergy"])
+rho = np.array(file["/PartType0/Densities"])
+P = np.array(file["/PartType0/Pressures"])
+u = np.array(file["/PartType0/InternalEnergies"])
m = np.array(file["/PartType0/Masses"])
vs = np.array(file["/PartType0/Velocities"])
ids = np.array(file["/PartType0/ParticleIDs"])
diff --git a/examples/HydroTests/SineWavePotential_2D/plotSolution.py b/examples/HydroTests/SineWavePotential_2D/plotSolution.py
index ee02f59c404db87a790465d2786e6296525e36b0..5c87b0f4f3682d486063522715517763b1035567 100644
--- a/examples/HydroTests/SineWavePotential_2D/plotSolution.py
+++ b/examples/HydroTests/SineWavePotential_2D/plotSolution.py
@@ -38,8 +38,8 @@ fileName = sys.argv[1]
file = h5py.File(fileName, 'r')
coords = np.array(file["/PartType0/Coordinates"])
-rho = np.array(file["/PartType0/Density"])
-u = np.array(file["/PartType0/InternalEnergy"])
+rho = np.array(file["/PartType0/Densities"])
+u = np.array(file["/PartType0/InternalEnergies"])
agrav = np.array(file["/PartType0/GravAcceleration"])
m = np.array(file["/PartType0/Masses"])
ids = np.array(file["/PartType0/ParticleIDs"])
diff --git a/examples/HydroTests/SineWavePotential_3D/plotSolution.py b/examples/HydroTests/SineWavePotential_3D/plotSolution.py
index 13cae037b64eff4ad4fec0003bf0f5ad3ce94896..7bfa82a5990572c478976614c50107e5254f0e00 100644
--- a/examples/HydroTests/SineWavePotential_3D/plotSolution.py
+++ b/examples/HydroTests/SineWavePotential_3D/plotSolution.py
@@ -38,8 +38,8 @@ fileName = sys.argv[1]
file = h5py.File(fileName, 'r')
coords = np.array(file["/PartType0/Coordinates"])
-rho = np.array(file["/PartType0/Density"])
-u = np.array(file["/PartType0/InternalEnergy"])
+rho = np.array(file["/PartType0/Densities"])
+u = np.array(file["/PartType0/InternalEnergies"])
agrav = np.array(file["/PartType0/GravAcceleration"])
m = np.array(file["/PartType0/Masses"])
ids = np.array(file["/PartType0/ParticleIDs"])
diff --git a/examples/HydroTests/SodShockSpherical_2D/plotSolution.py b/examples/HydroTests/SodShockSpherical_2D/plotSolution.py
index 57b7f7ddc64bc25df031eb0cba7547f40d46b31a..61060631eea1a7320ef207457d35031318bceccf 100644
--- a/examples/HydroTests/SodShockSpherical_2D/plotSolution.py
+++ b/examples/HydroTests/SodShockSpherical_2D/plotSolution.py
@@ -74,10 +74,10 @@ coords = sim["/PartType0/Coordinates"]
x = sqrt((coords[:,0] - 0.5)**2 + (coords[:,1] - 0.5)**2)
vels = sim["/PartType0/Velocities"]
v = sqrt(vels[:,0]**2 + vels[:,1]**2)
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
# Bin the data
rho_bin,x_bin_edge,_ = \
diff --git a/examples/HydroTests/SodShockSpherical_3D/plotSolution.py b/examples/HydroTests/SodShockSpherical_3D/plotSolution.py
index 539bfba799e3231bd26ae2eb39c271baa1fa6a4b..0a92f3aaf1831d67cb59dd71bc08cd8d973d9def 100644
--- a/examples/HydroTests/SodShockSpherical_3D/plotSolution.py
+++ b/examples/HydroTests/SodShockSpherical_3D/plotSolution.py
@@ -75,10 +75,10 @@ x = sqrt((coords[:,0] - 0.5)**2 + (coords[:,1] - 0.5)**2 + \
(coords[:,2] - 0.5)**2)
vels = sim["/PartType0/Velocities"]
v = sqrt(vels[:,0]**2 + vels[:,1]**2 + vels[:,2]**2)
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
# Bin the data
rho_bin,x_bin_edge,_ = \
diff --git a/examples/HydroTests/SodShock_1D/plotSolution.py b/examples/HydroTests/SodShock_1D/plotSolution.py
index a7e6d374bac616440dace666b85c3e7ade479bcd..770d05ac0493323c2cdd0f5905e409113d1a9eae 100644
--- a/examples/HydroTests/SodShock_1D/plotSolution.py
+++ b/examples/HydroTests/SodShock_1D/plotSolution.py
@@ -79,12 +79,12 @@ git = str(sim["Code"].attrs["Git Revision"])
x = sim["/PartType0/Coordinates"][:,0]
v = sim["/PartType0/Velocities"][:,0]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
try:
- alpha = sim["/PartType0/Viscosity"][:]
+ alpha = sim["/PartType0/ViscosityParameters"][:]
plot_alpha = True
except:
plot_alpha = False
diff --git a/examples/HydroTests/SodShock_2D/plotSolution.py b/examples/HydroTests/SodShock_2D/plotSolution.py
index 19cbe0ffb766845c051ffb6cea81bd918d890e36..769079da8824d58535e239bd8b54b592ce981a37 100644
--- a/examples/HydroTests/SodShock_2D/plotSolution.py
+++ b/examples/HydroTests/SodShock_2D/plotSolution.py
@@ -79,10 +79,10 @@ git = sim["Code"].attrs["Git Revision"]
x = sim["/PartType0/Coordinates"][:,0]
v = sim["/PartType0/Velocities"][:,0]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
N = 1000 # Number of points
x_min = -1.
diff --git a/examples/HydroTests/SodShock_3D/plotSolution.py b/examples/HydroTests/SodShock_3D/plotSolution.py
index 69b2fe4887e986156ed01e0f4177d01ccbed6035..c2028cedcea820e56c07962fe3ca2f1fe1347d40 100644
--- a/examples/HydroTests/SodShock_3D/plotSolution.py
+++ b/examples/HydroTests/SodShock_3D/plotSolution.py
@@ -79,19 +79,19 @@ git = sim["Code"].attrs["Git Revision"]
x = sim["/PartType0/Coordinates"][:,0]
v = sim["/PartType0/Velocities"][:,0]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
try:
- diffusion = sim["/PartType0/Diffusion"][:]
+ diffusion = sim["/PartType0/DiffusionParameters"][:]
plot_diffusion = True
except:
plot_diffusion = False
try:
- viscosity = sim["/PartType0/Viscosity"][:]
+ viscosity = sim["/PartType0/ViscosityParameters"][:]
plot_viscosity = True
except:
plot_viscosity = False
diff --git a/examples/HydroTests/SodShock_BCC_3D/plotSolution.py b/examples/HydroTests/SodShock_BCC_3D/plotSolution.py
index 365b679991e9a3a5bbb9e9d5108066c04e497c2f..9660e12a3d226bd6b0e3c152031c93cedb345933 100644
--- a/examples/HydroTests/SodShock_BCC_3D/plotSolution.py
+++ b/examples/HydroTests/SodShock_BCC_3D/plotSolution.py
@@ -94,10 +94,10 @@ time = sim.metadata.t.value
data = dict(
x=sim.gas.coordinates.value[:, 0],
v=sim.gas.velocities.value[:, 0],
- u=sim.gas.internal_energy.value,
- S=sim.gas.entropy.value,
- P=sim.gas.pressure.value,
- rho=sim.gas.density.value,
+ u=sim.gas.internal_energies.value,
+ S=sim.gas.entropies.value,
+ P=sim.gas.pressures.value,
+ rho=sim.gas.densities.value,
y=sim.gas.coordinates.value[:, 1],
z=sim.gas.coordinates.value[:, 2],
)
@@ -164,12 +164,9 @@ for key, label in plot.items():
zorder=-1,
)
- mask_noraster = np.logical_and.reduce([
- data["y"] < 0.52,
- data["y"] > 0.48,
- data["z"] < 0.52,
- data["z"] > 0.48
- ])
+ mask_noraster = np.logical_and.reduce(
+ [data["y"] < 0.52, data["y"] > 0.48, data["z"] < 0.52, data["z"] > 0.48]
+ )
axis.plot(
data["x"][mask_noraster],
diff --git a/examples/HydroTests/SquareTest_2D/plotSolutionLegacy.py b/examples/HydroTests/SquareTest_2D/plotSolutionLegacy.py
index 956da800c9096232d2e82cf4cff4c780672e0a8f..d8701c3d44390f1d2637f798c0e9af23531c4600 100644
--- a/examples/HydroTests/SquareTest_2D/plotSolutionLegacy.py
+++ b/examples/HydroTests/SquareTest_2D/plotSolutionLegacy.py
@@ -88,10 +88,10 @@ while centre_y < 0.:
pos = sim["/PartType0/Coordinates"][:,:]
vel = sim["/PartType0/Velocities"][:,:]
v_norm = sqrt(vel[:,0]**2 + vel[:,1]**2)
-rho = sim["/PartType0/Density"][:]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Densities"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
x = pos[:,0] - centre_x
y = pos[:,1] - centre_y
diff --git a/examples/HydroTests/VacuumSpherical_2D/plotSolution.py b/examples/HydroTests/VacuumSpherical_2D/plotSolution.py
index 6a65206ae20ccf79392054d047ba6be04f169f3e..de551c0ef4a5f4e027402c881069b7b4780f43d8 100644
--- a/examples/HydroTests/VacuumSpherical_2D/plotSolution.py
+++ b/examples/HydroTests/VacuumSpherical_2D/plotSolution.py
@@ -63,12 +63,12 @@ snap = int(sys.argv[1])
file = h5py.File("vacuum_{0:04d}.hdf5".format(snap), "r")
coords = file["/PartType0/Coordinates"]
x = np.sqrt((coords[:,0] - 0.5)**2 + (coords[:,1] - 0.5)**2)
-rho = file["/PartType0/Density"][:]
+rho = file["/PartType0/Densities"][:]
vels = file["/PartType0/Velocities"]
v = np.sqrt(vels[:,0]**2 + vels[:,1]**2)
-u = file["/PartType0/InternalEnergy"][:]
-S = file["/PartType0/Entropy"][:]
-P = file["/PartType0/Pressure"][:]
+u = file["/PartType0/InternalEnergies"][:]
+S = file["/PartType0/Entropies"][:]
+P = file["/PartType0/Pressures"][:]
time = file["/Header"].attrs["Time"][0]
scheme = file["/HydroScheme"].attrs["Scheme"]
diff --git a/examples/HydroTests/VacuumSpherical_3D/plotSolution.py b/examples/HydroTests/VacuumSpherical_3D/plotSolution.py
index c73e48ee2d311692cdf4aa3b0e52f4766b339df8..4a04acde575eae46de67c70b536b8774befd96ae 100644
--- a/examples/HydroTests/VacuumSpherical_3D/plotSolution.py
+++ b/examples/HydroTests/VacuumSpherical_3D/plotSolution.py
@@ -64,12 +64,12 @@ file = h5py.File("vacuum_{0:04d}.hdf5".format(snap), "r")
coords = file["/PartType0/Coordinates"]
x = np.sqrt((coords[:,0] - 0.5)**2 + (coords[:,1] - 0.5)**2 + \
(coords[:,2] - 0.5)**2)
-rho = file["/PartType0/Density"][:]
+rho = file["/PartType0/Densities"][:]
vels = file["/PartType0/Velocities"]
v = np.sqrt(vels[:,0]**2 + vels[:,1]**2 + vels[:,2]**2)
-u = file["/PartType0/InternalEnergy"][:]
-S = file["/PartType0/Entropy"][:]
-P = file["/PartType0/Pressure"][:]
+u = file["/PartType0/InternalEnergies"][:]
+S = file["/PartType0/Entropies"][:]
+P = file["/PartType0/Pressures"][:]
time = file["/Header"].attrs["Time"][0]
scheme = file["/HydroScheme"].attrs["Scheme"]
diff --git a/examples/HydroTests/Vacuum_1D/plotSolution.py b/examples/HydroTests/Vacuum_1D/plotSolution.py
index fceac10c25fd58b5bbcb6e31884cd62b4cfd61f5..eac7dc9e3ac43822ad167372f9f33bf2f5af0e2a 100644
--- a/examples/HydroTests/Vacuum_1D/plotSolution.py
+++ b/examples/HydroTests/Vacuum_1D/plotSolution.py
@@ -61,11 +61,11 @@ snap = int(sys.argv[1])
# Open the file and read the relevant data
file = h5py.File("vacuum_{0:04d}.hdf5".format(snap), "r")
x = file["/PartType0/Coordinates"][:,0]
-rho = file["/PartType0/Density"]
+rho = file["/PartType0/Densities"]
v = file["/PartType0/Velocities"][:,0]
-u = file["/PartType0/InternalEnergy"]
-S = file["/PartType0/Entropy"]
-P = file["/PartType0/Pressure"]
+u = file["/PartType0/InternalEnergies"]
+S = file["/PartType0/Entropies"]
+P = file["/PartType0/Pressures"]
time = file["/Header"].attrs["Time"][0]
scheme = file["/HydroScheme"].attrs["Scheme"]
diff --git a/examples/HydroTests/Vacuum_2D/plotSolution.py b/examples/HydroTests/Vacuum_2D/plotSolution.py
index 4d197234237df10b8cdbf197048a65991da023cf..ffd0eb1cdd857764f7ecc4e2d0c93fee3c5f29e8 100644
--- a/examples/HydroTests/Vacuum_2D/plotSolution.py
+++ b/examples/HydroTests/Vacuum_2D/plotSolution.py
@@ -62,11 +62,11 @@ snap = int(sys.argv[1])
# Open the file and read the relevant data
file = h5py.File("vacuum_{0:04d}.hdf5".format(snap), "r")
x = file["/PartType0/Coordinates"][:,0]
-rho = file["/PartType0/Density"][:]
+rho = file["/PartType0/Densities"][:]
v = file["/PartType0/Velocities"][:,0]
-u = file["/PartType0/InternalEnergy"][:]
-S = file["/PartType0/Entropy"][:]
-P = file["/PartType0/Pressure"][:]
+u = file["/PartType0/InternalEnergies"][:]
+S = file["/PartType0/Entropies"][:]
+P = file["/PartType0/Pressures"][:]
time = file["/Header"].attrs["Time"][0]
scheme = file["/HydroScheme"].attrs["Scheme"]
diff --git a/examples/HydroTests/Vacuum_3D/plotSolution.py b/examples/HydroTests/Vacuum_3D/plotSolution.py
index 4d197234237df10b8cdbf197048a65991da023cf..ffd0eb1cdd857764f7ecc4e2d0c93fee3c5f29e8 100644
--- a/examples/HydroTests/Vacuum_3D/plotSolution.py
+++ b/examples/HydroTests/Vacuum_3D/plotSolution.py
@@ -62,11 +62,11 @@ snap = int(sys.argv[1])
# Open the file and read the relevant data
file = h5py.File("vacuum_{0:04d}.hdf5".format(snap), "r")
x = file["/PartType0/Coordinates"][:,0]
-rho = file["/PartType0/Density"][:]
+rho = file["/PartType0/Densities"][:]
v = file["/PartType0/Velocities"][:,0]
-u = file["/PartType0/InternalEnergy"][:]
-S = file["/PartType0/Entropy"][:]
-P = file["/PartType0/Pressure"][:]
+u = file["/PartType0/InternalEnergies"][:]
+S = file["/PartType0/Entropies"][:]
+P = file["/PartType0/Pressures"][:]
time = file["/Header"].attrs["Time"][0]
scheme = file["/HydroScheme"].attrs["Scheme"]
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plotSolution.py b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plotSolution.py
index 89a87923148cb6872ab17b6d7229aef597ef3287..1ff8df3569f25590e5acb8046edeca0a1333d556 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plotSolution.py
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plotSolution.py
@@ -29,7 +29,7 @@ rcParams.update(params)
rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
snap = int(sys.argv[1])
-filename = "output_%.4d.hdf5"%snap
+filename = "output_%.4d.hdf5" % snap
f = h5.File(filename, "r")
@@ -40,7 +40,7 @@ year_in_cgs = 3600.0 * 24 * 365.0
Msun_in_cgs = 1.98848e33
G_in_cgs = 6.67259e-8
pc_in_cgs = 3.08567758e18
-Msun_p_pc2 = Msun_in_cgs / pc_in_cgs**2
+Msun_p_pc2 = Msun_in_cgs / pc_in_cgs ** 2
# Gemoetry info
boxsize = f["/Header"].attrs["BoxSize"]
@@ -52,66 +52,94 @@ unit_mass_in_cgs = f["/Units"].attrs["Unit mass in cgs (U_M)"]
unit_time_in_cgs = f["/Units"].attrs["Unit time in cgs (U_t)"]
# Calculate Gravitational constant in internal units
-G = G_in_cgs * ( unit_length_in_cgs**3 / unit_mass_in_cgs / unit_time_in_cgs**2)**(-1)
+G = G_in_cgs * (unit_length_in_cgs ** 3 / unit_mass_in_cgs / unit_time_in_cgs ** 2) ** (
+ -1
+)
# Read parameters of the SF model
KS_law_slope = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_exponent"])
KS_law_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_normalisation"])
KS_thresh_Z0 = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_Z0"])
KS_thresh_slope = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_slope"])
-KS_thresh_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_norm_H_p_cm3"])
+KS_thresh_norm = float(
+ f["/Parameters"].attrs["EAGLEStarFormation:threshold_norm_H_p_cm3"]
+)
KS_gas_fraction = float(f["/Parameters"].attrs["EAGLEStarFormation:gas_fraction"])
-KS_thresh_max_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_max_density_H_p_cm3"])
-KS_high_den_thresh = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_threshold_H_p_cm3"])
-KS_law_slope_high_den = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_exponent"])
-EOS_gamma_effective = float(f["/Parameters"].attrs["EAGLEStarFormation:EOS_gamma_effective"])
-EOS_density_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:EOS_density_norm_H_p_cm3"])
-EOS_temp_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:EOS_temperature_norm_K"])
+KS_thresh_max_norm = float(
+ f["/Parameters"].attrs["EAGLEStarFormation:threshold_max_density_H_p_cm3"]
+)
+KS_high_den_thresh = float(
+ f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_threshold_H_p_cm3"]
+)
+KS_law_slope_high_den = float(
+ f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_exponent"]
+)
+EOS_gamma_effective = float(
+ f["/Parameters"].attrs["EAGLEStarFormation:EOS_gamma_effective"]
+)
+EOS_density_norm = float(
+ f["/Parameters"].attrs["EAGLEStarFormation:EOS_density_norm_H_p_cm3"]
+)
+EOS_temp_norm = float(
+ f["/Parameters"].attrs["EAGLEStarFormation:EOS_temperature_norm_K"]
+)
# Read reference metallicity
EAGLE_Z = float(f["/Parameters"].attrs["EAGLEChemistry:init_abundance_metal"])
# Read parameters of the entropy floor
-EAGLEfloor_Jeans_rho_norm = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_density_threshold_H_p_cm3"])
-EAGLEfloor_Jeans_temperature_norm_K = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_temperature_norm_K"])
-EAGLEfloor_Jeans_gamma_effective = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_gamma_effective"])
-EAGLEfloor_cool_rho_norm = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_density_threshold_H_p_cm3"])
-EAGLEfloor_cool_temperature_norm_K = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_temperature_norm_K"])
-EAGLEfloor_cool_gamma_effective = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_gamma_effective"])
+EAGLEfloor_Jeans_rho_norm = float(
+ f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_density_threshold_H_p_cm3"]
+)
+EAGLEfloor_Jeans_temperature_norm_K = float(
+ f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_temperature_norm_K"]
+)
+EAGLEfloor_Jeans_gamma_effective = float(
+ f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_gamma_effective"]
+)
+EAGLEfloor_cool_rho_norm = float(
+ f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_density_threshold_H_p_cm3"]
+)
+EAGLEfloor_cool_temperature_norm_K = float(
+ f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_temperature_norm_K"]
+)
+EAGLEfloor_cool_gamma_effective = float(
+ f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_gamma_effective"]
+)
# Properties of the KS law
-KS_law_norm_cgs = KS_law_norm * Msun_in_cgs / ( 1e6 * pc_in_cgs**2 * year_in_cgs )
-gamma = 5./3.
+KS_law_norm_cgs = KS_law_norm * Msun_in_cgs / (1e6 * pc_in_cgs ** 2 * year_in_cgs)
+gamma = 5.0 / 3.0
EOS_press_norm = k_in_cgs * EOS_temp_norm * EOS_density_norm
# Star formation threshold
-SF_thresh = KS_thresh_norm * (EAGLE_Z / KS_thresh_Z0)**(KS_thresh_slope)
+SF_thresh = KS_thresh_norm * (EAGLE_Z / KS_thresh_Z0) ** (KS_thresh_slope)
# Read gas properties
gas_pos = f["/PartType0/Coordinates"][:, :]
gas_mass = f["/PartType0/Masses"][:]
-gas_rho = f["/PartType0/Density"][:]
+gas_rho = f["/PartType0/Densities"][:]
gas_T = f["/PartType0/Temperature"][:]
gas_SFR = f["/PartType0/SFR"][:]
-gas_XH = f["/PartType0/ElementAbundance"][:, 0]
-gas_Z = f["/PartType0/Metallicity"][:]
-gas_hsml = f["/PartType0/SmoothingLength"][:]
+gas_XH = f["/PartType0/ElementMassFractions"][:, 0]
+gas_Z = f["/PartType0/Metallicities"][:]
+gas_hsml = f["/PartType0/SmoothingLengths"][:]
gas_sSFR = gas_SFR / gas_mass
# Read the Star properties
stars_pos = f["/PartType4/Coordinates"][:, :]
stars_BirthDensity = f["/PartType4/BirthDensity"][:]
stars_BirthTime = f["/PartType4/BirthTime"][:]
-stars_XH = f["/PartType4/ElementAbundance"][:,0]
+stars_XH = f["/PartType4/ElementAbundance"][:, 0]
# Centre the box
gas_pos[:, 0] -= centre[0]
gas_pos[:, 1] -= centre[1]
gas_pos[:, 2] -= centre[2]
-stars_pos[:,0] -= centre[0]
-stars_pos[:,1] -= centre[1]
-stars_pos[:,2] -= centre[2]
+stars_pos[:, 0] -= centre[0]
+stars_pos[:, 1] -= centre[1]
+stars_pos[:, 2] -= centre[2]
# Turn the mass into better units
gas_mass *= unit_mass_in_cgs / Msun_in_cgs
@@ -132,9 +160,13 @@ stars_BirthDensity *= stars_XH
# Equations of state
eos_cool_rho = np.logspace(-5, 5, 1000)
-eos_cool_T = EAGLEfloor_cool_temperature_norm_K * (eos_cool_rho / EAGLEfloor_cool_rho_norm) ** ( EAGLEfloor_cool_gamma_effective - 1.0 )
+eos_cool_T = EAGLEfloor_cool_temperature_norm_K * (
+ eos_cool_rho / EAGLEfloor_cool_rho_norm
+) ** (EAGLEfloor_cool_gamma_effective - 1.0)
eos_Jeans_rho = np.logspace(-1, 5, 1000)
-eos_Jeans_T = EAGLEfloor_Jeans_temperature_norm_K * (eos_Jeans_rho / EAGLEfloor_Jeans_rho_norm) ** (EAGLEfloor_Jeans_gamma_effective - 1.0 )
+eos_Jeans_T = EAGLEfloor_Jeans_temperature_norm_K * (
+ eos_Jeans_rho / EAGLEfloor_Jeans_rho_norm
+) ** (EAGLEfloor_Jeans_gamma_effective - 1.0)
########################################################################3
@@ -156,7 +188,15 @@ subplot(111, xscale="log", yscale="log")
plot(eos_cool_rho, eos_cool_T, "k--", lw=0.6)
plot(eos_Jeans_rho, eos_Jeans_T, "k--", lw=0.6)
plot([SF_thresh, SF_thresh], [1, 1e10], "k:", lw=0.6)
-text(SF_thresh*0.9, 2e4, "$n_{\\rm H, thresh}=%.3f~{\\rm cm^{-3}}$"%SF_thresh, fontsize=8, rotation=90, ha="right", va="bottom")
+text(
+ SF_thresh * 0.9,
+ 2e4,
+ "$n_{\\rm H, thresh}=%.3f~{\\rm cm^{-3}}$" % SF_thresh,
+ fontsize=8,
+ rotation=90,
+ ha="right",
+ va="bottom",
+)
scatter(gas_nH[gas_SFR > 0.0], gas_T[gas_SFR > 0.0], s=0.2)
xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=2)
@@ -188,37 +228,61 @@ star_mask = (
& (stars_pos[:, 2] > -1.0)
)
-stars_BirthDensity = stars_BirthDensity[star_mask]
-#stars_BirthFlag = stars_BirthFlag[star_mask]
+stars_BirthDensity = stars_BirthDensity[star_mask]
+# stars_BirthFlag = stars_BirthFlag[star_mask]
stars_BirthTime = stars_BirthTime[star_mask]
# Histogram of the birth density
figure()
subplot(111, xscale="linear", yscale="linear")
-hist(np.log10(stars_BirthDensity),density=True,bins=20,range=[-2,5])
+hist(np.log10(stars_BirthDensity), density=True, bins=20, range=[-2, 5])
xlabel("${\\rm Stellar~birth~density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
ylabel("${\\rm Probability}$", labelpad=-7)
savefig("BirthDensity.png", dpi=200)
# Plot of the specific star formation rate in the galaxy
-rhos = 10**np.linspace(-1,np.log10(KS_high_den_thresh),100)
-rhoshigh = 10**np.linspace(np.log10(KS_high_den_thresh),5,100)
+rhos = 10 ** np.linspace(-1, np.log10(KS_high_den_thresh), 100)
+rhoshigh = 10 ** np.linspace(np.log10(KS_high_den_thresh), 5, 100)
-P_effective = EOS_press_norm * ( rhos / EOS_density_norm)**(EOS_gamma_effective)
-P_norm_high = EOS_press_norm * (KS_high_den_thresh / EOS_density_norm)**(EOS_gamma_effective)
-sSFR = KS_law_norm_cgs * (Msun_p_pc2)**(-KS_law_slope) * (gamma/G_in_cgs * KS_gas_fraction *P_effective)**((KS_law_slope-1.)/2.)
-KS_law_norm_high_den_cgs = KS_law_norm_cgs * (Msun_p_pc2)**(-KS_law_slope) * (gamma/G_in_cgs * KS_gas_fraction * P_norm_high)**((KS_law_slope-1.)/2.)
-sSFR_high_den = KS_law_norm_high_den_cgs * ((rhoshigh/KS_high_den_thresh)**EOS_gamma_effective)**((KS_law_slope_high_den-1)/2.)
+P_effective = EOS_press_norm * (rhos / EOS_density_norm) ** (EOS_gamma_effective)
+P_norm_high = EOS_press_norm * (KS_high_den_thresh / EOS_density_norm) ** (
+ EOS_gamma_effective
+)
+sSFR = (
+ KS_law_norm_cgs
+ * (Msun_p_pc2) ** (-KS_law_slope)
+ * (gamma / G_in_cgs * KS_gas_fraction * P_effective) ** ((KS_law_slope - 1.0) / 2.0)
+)
+KS_law_norm_high_den_cgs = (
+ KS_law_norm_cgs
+ * (Msun_p_pc2) ** (-KS_law_slope)
+ * (gamma / G_in_cgs * KS_gas_fraction * P_norm_high) ** ((KS_law_slope - 1.0) / 2.0)
+)
+sSFR_high_den = KS_law_norm_high_den_cgs * (
+ (rhoshigh / KS_high_den_thresh) ** EOS_gamma_effective
+) ** ((KS_law_slope_high_den - 1) / 2.0)
# density - sSFR plane
figure()
subplot(111)
-hist2d(np.log10(gas_nH), np.log10(gas_sSFR), bins=50,range=[[-1.5,5],[-.5,2.5]])
-plot(np.log10(rhos),np.log10(sSFR)+np.log10(year_in_cgs)+9.,'k--',label='sSFR low density EAGLE')
-plot(np.log10(rhoshigh),np.log10(sSFR_high_den)+np.log10(year_in_cgs)+9.,'k--',label='sSFR high density EAGLE')
+hist2d(np.log10(gas_nH), np.log10(gas_sSFR), bins=50, range=[[-1.5, 5], [-0.5, 2.5]])
+plot(
+ np.log10(rhos),
+ np.log10(sSFR) + np.log10(year_in_cgs) + 9.0,
+ "k--",
+ label="sSFR low density EAGLE",
+)
+plot(
+ np.log10(rhoshigh),
+ np.log10(sSFR_high_den) + np.log10(year_in_cgs) + 9.0,
+ "k--",
+ label="sSFR high density EAGLE",
+)
xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=2)
ylabel("${\\rm sSFR}~[{\\rm Gyr^{-1}}]$", labelpad=0)
-xticks([-1, 0, 1, 2, 3, 4], ["$10^{-1}$", "$10^0$", "$10^1$", "$10^2$", "$10^3$", "$10^4$"])
+xticks(
+ [-1, 0, 1, 2, 3, 4], ["$10^{-1}$", "$10^0$", "$10^1$", "$10^2$", "$10^3$", "$10^4$"]
+)
yticks([0, 1, 2], ["$10^0$", "$10^1$", "$10^2$"])
xlim(-1.4, 4.9)
ylim(-0.5, 2.2)
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plot_box_evolution.py b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plot_box_evolution.py
index 67da3c390be1240323941b835e056dcd6e27feed..94f27c87ff46c48ffc6b0df8c3e02c7abb6df875 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plot_box_evolution.py
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plot_box_evolution.py
@@ -1,24 +1,25 @@
###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
- # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- #
- # 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/>.
- #
- ##############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+# Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#
+# 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/>.
+#
+##############################################################################
import matplotlib
+
matplotlib.use("Agg")
from pylab import *
from scipy import stats
@@ -27,40 +28,42 @@ import numpy as np
import glob
import os.path
-def find_indices(a,b):
- result = np.zeros(len(b))
- for i in range(len(b)):
- result[i] = ((np.where(a == b[i]))[0])[0]
- return result
+def find_indices(a, b):
+ result = np.zeros(len(b))
+ for i in range(len(b)):
+ result[i] = ((np.where(a == b[i]))[0])[0]
+
+ return result
# Plot parameters
-params = {'axes.labelsize': 10,
-'axes.titlesize': 10,
-'font.size': 12,
-'legend.fontsize': 12,
-'xtick.labelsize': 10,
-'ytick.labelsize': 10,
-'text.usetex': True,
- 'figure.figsize' : (9.90,6.45),
-'figure.subplot.left' : 0.1,
-'figure.subplot.right' : 0.99,
-'figure.subplot.bottom' : 0.1,
-'figure.subplot.top' : 0.95,
-'figure.subplot.wspace' : 0.2,
-'figure.subplot.hspace' : 0.2,
-'lines.markersize' : 6,
-'lines.linewidth' : 3.,
-'text.latex.unicode': True
+params = {
+ "axes.labelsize": 10,
+ "axes.titlesize": 10,
+ "font.size": 12,
+ "legend.fontsize": 12,
+ "xtick.labelsize": 10,
+ "ytick.labelsize": 10,
+ "text.usetex": True,
+ "figure.figsize": (9.90, 6.45),
+ "figure.subplot.left": 0.1,
+ "figure.subplot.right": 0.99,
+ "figure.subplot.bottom": 0.1,
+ "figure.subplot.top": 0.95,
+ "figure.subplot.wspace": 0.2,
+ "figure.subplot.hspace": 0.2,
+ "lines.markersize": 6,
+ "lines.linewidth": 3.0,
+ "text.latex.unicode": True,
}
rcParams.update(params)
-rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
# Number of snapshots and elements
-newest_snap_name = max(glob.glob('output_*.hdf5'), key=os.path.getctime)
-n_snapshots = int(newest_snap_name.replace('output_','').replace('.hdf5','')) + 1
+newest_snap_name = max(glob.glob("output_*.hdf5"), key=os.path.getctime)
+n_snapshots = int(newest_snap_name.replace("output_", "").replace(".hdf5", "")) + 1
n_elements = 9
# Read the simulation data
@@ -84,10 +87,10 @@ unit_energy_in_cgs = unit_mass_in_cgs * unit_vel_in_cgs * unit_vel_in_cgs
unit_length_in_si = 0.01 * unit_length_in_cgs
unit_mass_in_si = 0.001 * unit_mass_in_cgs
unit_time_in_si = unit_time_in_cgs
-unit_density_in_cgs = unit_mass_in_cgs*unit_length_in_cgs**-3
-unit_pressure_in_cgs = unit_mass_in_cgs/unit_length_in_cgs*unit_time_in_cgs**-2
-unit_int_energy_in_cgs = unit_energy_in_cgs/unit_mass_in_cgs
-unit_entropy_in_cgs = unit_energy_in_cgs/unit_temp_in_cgs
+unit_density_in_cgs = unit_mass_in_cgs * unit_length_in_cgs ** -3
+unit_pressure_in_cgs = unit_mass_in_cgs / unit_length_in_cgs * unit_time_in_cgs ** -2
+unit_int_energy_in_cgs = unit_energy_in_cgs / unit_mass_in_cgs
+unit_entropy_in_cgs = unit_energy_in_cgs / unit_temp_in_cgs
Gyr_in_cgs = 3.155e16
Msun_in_cgs = 1.989e33
@@ -95,61 +98,111 @@ box_energy = zeros(n_snapshots)
box_mass = zeros(n_snapshots)
box_star_mass = zeros(n_snapshots)
box_metal_mass = zeros(n_snapshots)
-element_mass = zeros((n_snapshots,n_elements))
+element_mass = zeros((n_snapshots, n_elements))
t = zeros(n_snapshots)
# Read data from snapshots
for i in range(n_snapshots):
- print("reading snapshot "+str(i))
- # Read the simulation data
- sim = h5py.File("output_%04d.hdf5"%i, "r")
- t[i] = sim["/Header"].attrs["Time"][0]
- #ids = sim["/PartType0/ParticleIDs"][:]
-
- masses = sim["/PartType0/Masses"][:]
- box_mass[i] = np.sum(masses)
-
- star_masses = sim["/PartType4/Masses"][:]
- box_star_mass[i] = np.sum(star_masses)
-
- metallicities = sim["/PartType0/Metallicity"][:]
- box_metal_mass[i] = np.sum(metallicities * masses)
-
- internal_energies = sim["/PartType0/InternalEnergy"][:]
- box_energy[i] = np.sum(masses * internal_energies)
+ print("reading snapshot " + str(i))
+ # Read the simulation data
+ sim = h5py.File("output_%04d.hdf5" % i, "r")
+ t[i] = sim["/Header"].attrs["Time"][0]
+ # ids = sim["/PartType0/ParticleIDs"][:]
+
+ masses = sim["/PartType0/Masses"][:]
+ box_mass[i] = np.sum(masses)
+
+ star_masses = sim["/PartType4/Masses"][:]
+ box_star_mass[i] = np.sum(star_masses)
+
+ metallicities = sim["/PartType0/Metallicities"][:]
+ box_metal_mass[i] = np.sum(metallicities * masses)
+
+ internal_energies = sim["/PartType0/InternalEnergies"][:]
+ box_energy[i] = np.sum(masses * internal_energies)
# Plot the interesting quantities
figure()
# Box mass --------------------------------
subplot(221)
-plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_mass[1:] - box_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', marker = "*", ms=0.5, label='swift')
+plot(
+ t[1:] * unit_time_in_cgs / Gyr_in_cgs,
+ (box_mass[1:] - box_mass[0]) * unit_mass_in_cgs / Msun_in_cgs,
+ linewidth=0.5,
+ color="k",
+ marker="*",
+ ms=0.5,
+ label="swift",
+)
xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
ylabel("Change in total gas particle mass (Msun)", labelpad=2)
-ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
# Box metal mass --------------------------------
subplot(222)
-plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_metal_mass[1:] - box_metal_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', ms=0.5, label='swift')
+plot(
+ t[1:] * unit_time_in_cgs / Gyr_in_cgs,
+ (box_metal_mass[1:] - box_metal_mass[0]) * unit_mass_in_cgs / Msun_in_cgs,
+ linewidth=0.5,
+ color="k",
+ ms=0.5,
+ label="swift",
+)
xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
ylabel("Change in total metal mass of gas particles (Msun)", labelpad=2)
-ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
# Box energy --------------------------------
subplot(223)
-plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_energy[1:] - box_energy[0])* unit_energy_in_cgs, linewidth=0.5, color='k', ms=0.5, label='swift')
+plot(
+ t[1:] * unit_time_in_cgs / Gyr_in_cgs,
+ (box_energy[1:] - box_energy[0]) * unit_energy_in_cgs,
+ linewidth=0.5,
+ color="k",
+ ms=0.5,
+ label="swift",
+)
xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
ylabel("Change in total energy of gas particles (erg)", labelpad=2)
-ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
# Box mass --------------------------------
subplot(224)
-plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_mass[1:] - box_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', marker = "*", ms=0.5, label='gas')
-plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_star_mass[1:] - box_star_mass[n_snapshots-1])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='r', marker = "*", ms=0.5, label='stars')
-plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_star_mass[1:] - box_star_mass[n_snapshots-1] + box_mass[1:] - box_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='g', marker = "*", ms=0.5, label='total')
+plot(
+ t[1:] * unit_time_in_cgs / Gyr_in_cgs,
+ (box_mass[1:] - box_mass[0]) * unit_mass_in_cgs / Msun_in_cgs,
+ linewidth=0.5,
+ color="k",
+ marker="*",
+ ms=0.5,
+ label="gas",
+)
+plot(
+ t[1:] * unit_time_in_cgs / Gyr_in_cgs,
+ (box_star_mass[1:] - box_star_mass[n_snapshots - 1])
+ * unit_mass_in_cgs
+ / Msun_in_cgs,
+ linewidth=0.5,
+ color="r",
+ marker="*",
+ ms=0.5,
+ label="stars",
+)
+plot(
+ t[1:] * unit_time_in_cgs / Gyr_in_cgs,
+ (box_star_mass[1:] - box_star_mass[n_snapshots - 1] + box_mass[1:] - box_mass[0])
+ * unit_mass_in_cgs
+ / Msun_in_cgs,
+ linewidth=0.5,
+ color="g",
+ marker="*",
+ ms=0.5,
+ label="total",
+)
xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
ylabel("Change in total gas particle mass (Msun)", labelpad=2)
-ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
legend()
savefig("box_evolution.png", dpi=200)
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/plotSolution.py b/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/plotSolution.py
index 73e4878e8e00a35fe19c359652be0d57153dea62..044ad2bc78958cbf669c7257122b1ff80a94ba1a 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/plotSolution.py
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/plotSolution.py
@@ -49,7 +49,7 @@ rcParams.update(params)
rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
snap = int(sys.argv[1])
-filename = "output_%.4d.hdf5"%snap
+filename = "output_%.4d.hdf5" % snap
f = h5.File(filename, "r")
@@ -60,7 +60,7 @@ year_in_cgs = 3600.0 * 24 * 365.0
Msun_in_cgs = 1.98848e33
G_in_cgs = 6.67259e-8
pc_in_cgs = 3.08567758e18
-Msun_p_pc2 = Msun_in_cgs / pc_in_cgs**2
+Msun_p_pc2 = Msun_in_cgs / pc_in_cgs ** 2
# Gemoetry info
boxsize = f["/Header"].attrs["BoxSize"]
@@ -72,57 +72,85 @@ unit_mass_in_cgs = f["/Units"].attrs["Unit mass in cgs (U_M)"]
unit_time_in_cgs = f["/Units"].attrs["Unit time in cgs (U_t)"]
# Calculate Gravitational constant in internal units
-G = G_in_cgs * ( unit_length_in_cgs**3 / unit_mass_in_cgs / unit_time_in_cgs**2)**(-1)
+G = G_in_cgs * (unit_length_in_cgs ** 3 / unit_mass_in_cgs / unit_time_in_cgs ** 2) ** (
+ -1
+)
# Read parameters of the SF model
KS_law_slope = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_exponent"])
KS_law_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_normalisation"])
KS_thresh_Z0 = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_Z0"])
KS_thresh_slope = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_slope"])
-KS_thresh_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_norm_H_p_cm3"])
+KS_thresh_norm = float(
+ f["/Parameters"].attrs["EAGLEStarFormation:threshold_norm_H_p_cm3"]
+)
KS_gas_fraction = float(f["/Parameters"].attrs["EAGLEStarFormation:gas_fraction"])
-KS_thresh_max_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_max_density_H_p_cm3"])
-KS_high_den_thresh = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_threshold_H_p_cm3"])
-KS_law_slope_high_den = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_exponent"])
-EOS_gamma_effective = float(f["/Parameters"].attrs["EAGLEStarFormation:EOS_gamma_effective"])
-EOS_density_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:EOS_density_norm_H_p_cm3"])
-EOS_temp_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:EOS_temperature_norm_K"])
+KS_thresh_max_norm = float(
+ f["/Parameters"].attrs["EAGLEStarFormation:threshold_max_density_H_p_cm3"]
+)
+KS_high_den_thresh = float(
+ f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_threshold_H_p_cm3"]
+)
+KS_law_slope_high_den = float(
+ f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_exponent"]
+)
+EOS_gamma_effective = float(
+ f["/Parameters"].attrs["EAGLEStarFormation:EOS_gamma_effective"]
+)
+EOS_density_norm = float(
+ f["/Parameters"].attrs["EAGLEStarFormation:EOS_density_norm_H_p_cm3"]
+)
+EOS_temp_norm = float(
+ f["/Parameters"].attrs["EAGLEStarFormation:EOS_temperature_norm_K"]
+)
# Read reference metallicity
EAGLE_Z = float(f["/Parameters"].attrs["EAGLEChemistry:init_abundance_metal"])
# Read parameters of the entropy floor
-EAGLEfloor_Jeans_rho_norm = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_density_threshold_H_p_cm3"])
-EAGLEfloor_Jeans_temperature_norm_K = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_temperature_norm_K"])
-EAGLEfloor_Jeans_gamma_effective = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_gamma_effective"])
-EAGLEfloor_cool_rho_norm = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_density_threshold_H_p_cm3"])
-EAGLEfloor_cool_temperature_norm_K = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_temperature_norm_K"])
-EAGLEfloor_cool_gamma_effective = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_gamma_effective"])
+EAGLEfloor_Jeans_rho_norm = float(
+ f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_density_threshold_H_p_cm3"]
+)
+EAGLEfloor_Jeans_temperature_norm_K = float(
+ f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_temperature_norm_K"]
+)
+EAGLEfloor_Jeans_gamma_effective = float(
+ f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_gamma_effective"]
+)
+EAGLEfloor_cool_rho_norm = float(
+ f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_density_threshold_H_p_cm3"]
+)
+EAGLEfloor_cool_temperature_norm_K = float(
+ f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_temperature_norm_K"]
+)
+EAGLEfloor_cool_gamma_effective = float(
+ f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_gamma_effective"]
+)
# Properties of the KS law
-KS_law_norm_cgs = KS_law_norm * Msun_in_cgs / ( 1e6 * pc_in_cgs**2 * year_in_cgs )
-gamma = 5./3.
+KS_law_norm_cgs = KS_law_norm * Msun_in_cgs / (1e6 * pc_in_cgs ** 2 * year_in_cgs)
+gamma = 5.0 / 3.0
EOS_press_norm = k_in_cgs * EOS_temp_norm * EOS_density_norm
# Star formation threshold
-SF_thresh = KS_thresh_norm * (EAGLE_Z / KS_thresh_Z0)**(KS_thresh_slope)
+SF_thresh = KS_thresh_norm * (EAGLE_Z / KS_thresh_Z0) ** (KS_thresh_slope)
# Read gas properties
gas_pos = f["/PartType0/Coordinates"][:, :]
gas_mass = f["/PartType0/Masses"][:]
-gas_rho = f["/PartType0/Density"][:]
+gas_rho = f["/PartType0/Densities"][:]
gas_T = f["/PartType0/Temperature"][:]
gas_SFR = f["/PartType0/SFR"][:]
-gas_XH = f["/PartType0/ElementAbundance"][:, 0]
-gas_Z = f["/PartType0/Metallicity"][:]
-gas_hsml = f["/PartType0/SmoothingLength"][:]
+gas_XH = f["/PartType0/ElementMassFractions"][:, 0]
+gas_Z = f["/PartType0/Metallicities"][:]
+gas_hsml = f["/PartType0/SmoothingLengths"][:]
gas_sSFR = gas_SFR / gas_mass
# Read the Star properties
stars_pos = f["/PartType4/Coordinates"][:, :]
stars_BirthDensity = f["/PartType4/BirthDensity"][:]
stars_BirthTime = f["/PartType4/BirthTime"][:]
-stars_XH = f["/PartType4/ElementAbundance"][:,0]
+stars_XH = f["/PartType4/ElementAbundance"][:, 0]
stars_BirthTemperature = f["/PartType4/BirthTemperature"][:]
# Centre the box
@@ -130,9 +158,9 @@ gas_pos[:, 0] -= centre[0]
gas_pos[:, 1] -= centre[1]
gas_pos[:, 2] -= centre[2]
-stars_pos[:,0] -= centre[0]
-stars_pos[:,1] -= centre[1]
-stars_pos[:,2] -= centre[2]
+stars_pos[:, 0] -= centre[0]
+stars_pos[:, 1] -= centre[1]
+stars_pos[:, 2] -= centre[2]
# Turn the mass into better units
gas_mass *= unit_mass_in_cgs / Msun_in_cgs
@@ -156,9 +184,13 @@ stars_BirthDensity *= stars_XH
# Equations of state
eos_cool_rho = np.logspace(-5, 5, 1000)
-eos_cool_T = EAGLEfloor_cool_temperature_norm_K * (eos_cool_rho / EAGLEfloor_cool_rho_norm) ** ( EAGLEfloor_cool_gamma_effective - 1.0 )
+eos_cool_T = EAGLEfloor_cool_temperature_norm_K * (
+ eos_cool_rho / EAGLEfloor_cool_rho_norm
+) ** (EAGLEfloor_cool_gamma_effective - 1.0)
eos_Jeans_rho = np.logspace(-1, 5, 1000)
-eos_Jeans_T = EAGLEfloor_Jeans_temperature_norm_K * (eos_Jeans_rho / EAGLEfloor_Jeans_rho_norm) ** (EAGLEfloor_Jeans_gamma_effective - 1.0 )
+eos_Jeans_T = EAGLEfloor_Jeans_temperature_norm_K * (
+ eos_Jeans_rho / EAGLEfloor_Jeans_rho_norm
+) ** (EAGLEfloor_Jeans_gamma_effective - 1.0)
########################################################################3
@@ -180,7 +212,15 @@ subplot(111, xscale="log", yscale="log")
plot(eos_cool_rho, eos_cool_T, "k--", lw=0.6)
plot(eos_Jeans_rho, eos_Jeans_T, "k--", lw=0.6)
plot([SF_thresh, SF_thresh], [1, 1e10], "k:", lw=0.6)
-text(SF_thresh*0.9, 2e4, "$n_{\\rm H, thresh}=%.3f~{\\rm cm^{-3}}$"%SF_thresh, fontsize=8, rotation=90, ha="right", va="bottom")
+text(
+ SF_thresh * 0.9,
+ 2e4,
+ "$n_{\\rm H, thresh}=%.3f~{\\rm cm^{-3}}$" % SF_thresh,
+ fontsize=8,
+ rotation=90,
+ ha="right",
+ va="bottom",
+)
scatter(gas_nH[gas_SFR > 0.0], gas_T[gas_SFR > 0.0], s=0.2)
xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=2)
@@ -199,66 +239,92 @@ star_mask = (
& (stars_pos[:, 2] > -1.0)
)
-stars_BirthDensity = stars_BirthDensity[star_mask]
-#stars_BirthFlag = stars_BirthFlag[star_mask]
+stars_BirthDensity = stars_BirthDensity[star_mask]
+# stars_BirthFlag = stars_BirthFlag[star_mask]
stars_BirthTime = stars_BirthTime[star_mask]
# Histogram of the birth density
figure()
subplot(111, xscale="linear", yscale="linear")
-hist(np.log10(stars_BirthDensity),density=True,bins=20,range=[-2,5])
+hist(np.log10(stars_BirthDensity), density=True, bins=20, range=[-2, 5])
xlabel("${\\rm Stellar~birth~density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
ylabel("${\\rm Probability}$", labelpad=3)
savefig("BirthDensity.png", dpi=200)
-# Histogram of the birth temperature
+# Histogram of the birth temperature
figure()
subplot(111, xscale="linear", yscale="linear")
-hist(np.log10(stars_BirthTemperature),density=True,bins=20,range=[3.5,5.0])
+hist(np.log10(stars_BirthTemperature), density=True, bins=20, range=[3.5, 5.0])
xlabel("${\\rm Stellar~birth~temperature}~[{\\rm K}]$", labelpad=0)
ylabel("${\\rm Probability}$", labelpad=3)
savefig("BirthTemperature.png", dpi=200)
# Plot of the specific star formation rate in the galaxy
-rhos = 10**np.linspace(-1,np.log10(KS_high_den_thresh),100)
-rhoshigh = 10**np.linspace(np.log10(KS_high_den_thresh),5,100)
+rhos = 10 ** np.linspace(-1, np.log10(KS_high_den_thresh), 100)
+rhoshigh = 10 ** np.linspace(np.log10(KS_high_den_thresh), 5, 100)
-P_effective = EOS_press_norm * ( rhos / EOS_density_norm)**(EOS_gamma_effective)
-P_norm_high = EOS_press_norm * (KS_high_den_thresh / EOS_density_norm)**(EOS_gamma_effective)
-sSFR = KS_law_norm_cgs * (Msun_p_pc2)**(-KS_law_slope) * (gamma/G_in_cgs * KS_gas_fraction *P_effective)**((KS_law_slope-1.)/2.)
-KS_law_norm_high_den_cgs = KS_law_norm_cgs * (Msun_p_pc2)**(-KS_law_slope) * (gamma/G_in_cgs * KS_gas_fraction * P_norm_high)**((KS_law_slope-1.)/2.)
-sSFR_high_den = KS_law_norm_high_den_cgs * ((rhoshigh/KS_high_den_thresh)**EOS_gamma_effective)**((KS_law_slope_high_den-1)/2.)
+P_effective = EOS_press_norm * (rhos / EOS_density_norm) ** (EOS_gamma_effective)
+P_norm_high = EOS_press_norm * (KS_high_den_thresh / EOS_density_norm) ** (
+ EOS_gamma_effective
+)
+sSFR = (
+ KS_law_norm_cgs
+ * (Msun_p_pc2) ** (-KS_law_slope)
+ * (gamma / G_in_cgs * KS_gas_fraction * P_effective) ** ((KS_law_slope - 1.0) / 2.0)
+)
+KS_law_norm_high_den_cgs = (
+ KS_law_norm_cgs
+ * (Msun_p_pc2) ** (-KS_law_slope)
+ * (gamma / G_in_cgs * KS_gas_fraction * P_norm_high) ** ((KS_law_slope - 1.0) / 2.0)
+)
+sSFR_high_den = KS_law_norm_high_den_cgs * (
+ (rhoshigh / KS_high_den_thresh) ** EOS_gamma_effective
+) ** ((KS_law_slope_high_den - 1) / 2.0)
# density - sSFR plane
figure()
subplot(111)
-hist2d(np.log10(gas_nH), np.log10(gas_sSFR), bins=50,range=[[-1.5,5],[-.5,2.5]])
-plot(np.log10(rhos),np.log10(sSFR)+np.log10(year_in_cgs)+9.,'k--',label='sSFR low density EAGLE')
-plot(np.log10(rhoshigh),np.log10(sSFR_high_den)+np.log10(year_in_cgs)+9.,'k--',label='sSFR high density EAGLE')
+hist2d(np.log10(gas_nH), np.log10(gas_sSFR), bins=50, range=[[-1.5, 5], [-0.5, 2.5]])
+plot(
+ np.log10(rhos),
+ np.log10(sSFR) + np.log10(year_in_cgs) + 9.0,
+ "k--",
+ label="sSFR low density EAGLE",
+)
+plot(
+ np.log10(rhoshigh),
+ np.log10(sSFR_high_den) + np.log10(year_in_cgs) + 9.0,
+ "k--",
+ label="sSFR high density EAGLE",
+)
xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=2)
ylabel("${\\rm sSFR}~[{\\rm Gyr^{-1}}]$", labelpad=0)
-xticks([-1, 0, 1, 2, 3, 4], ["$10^{-1}$", "$10^0$", "$10^1$", "$10^2$", "$10^3$", "$10^4$"])
+xticks(
+ [-1, 0, 1, 2, 3, 4], ["$10^{-1}$", "$10^0$", "$10^1$", "$10^2$", "$10^3$", "$10^4$"]
+)
yticks([0, 1, 2], ["$10^0$", "$10^1$", "$10^2$"])
xlim(-1.4, 4.9)
ylim(-0.5, 2.2)
savefig("density-sSFR.png", dpi=200)
-SFR_low = 10**(np.log10(sSFR)+np.log10(year_in_cgs)+np.log10(median_gas_mass))
-SFR_high = 10**(np.log10(sSFR_high_den)+np.log10(year_in_cgs)+np.log10(median_gas_mass))
-SFR_low_min = np.floor(np.log10(.75*np.min(SFR_low)))
-SFR_high_max = np.ceil(np.log10(1.25*np.max(SFR_high)))
+SFR_low = 10 ** (np.log10(sSFR) + np.log10(year_in_cgs) + np.log10(median_gas_mass))
+SFR_high = 10 ** (
+ np.log10(sSFR_high_den) + np.log10(year_in_cgs) + np.log10(median_gas_mass)
+)
+SFR_low_min = np.floor(np.log10(0.75 * np.min(SFR_low)))
+SFR_high_max = np.ceil(np.log10(1.25 * np.max(SFR_high)))
# 3D Density vs SFR
rcParams.update({"figure.subplot.left": 0.18})
figure()
subplot(111, xscale="log", yscale="log")
scatter(gas_nH, gas_SFR, s=0.2)
-plot(rhos,SFR_low,'k--',lw=1,label='SFR low density EAGLE')
-plot(rhoshigh,SFR_high,'k--',lw=1,label='SFR high density EAGLE')
+plot(rhos, SFR_low, "k--", lw=1, label="SFR low density EAGLE")
+plot(rhoshigh, SFR_high, "k--", lw=1, label="SFR high density EAGLE")
xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
ylabel("${\\rm SFR}~[{\\rm M_\\odot~\\cdot~yr^{-1}}]$", labelpad=2)
xlim(1e-2, 1e5)
-ylim(10**SFR_low_min, 10**(SFR_high_max+0.1))
+ylim(10 ** SFR_low_min, 10 ** (SFR_high_max + 0.1))
savefig("rho_SFR.png", dpi=200)
rcParams.update({"figure.subplot.left": 0.15})
########################################################################3
diff --git a/examples/SantaBarbara/SantaBarbara-256/plotTempEvolution.py b/examples/SantaBarbara/SantaBarbara-256/plotTempEvolution.py
index dab4b2c90a7b751c8d143ed38c614473c951988a..63e46ccaee7be9ea18090e13ae15bb0a1fae4bef 100644
--- a/examples/SantaBarbara/SantaBarbara-256/plotTempEvolution.py
+++ b/examples/SantaBarbara/SantaBarbara-256/plotTempEvolution.py
@@ -128,7 +128,7 @@ for i in range(n_snapshots):
z[i] = sim["/Cosmology"].attrs["Redshift"][0]
a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
- u = sim["/PartType0/InternalEnergy"][:]
+ u = sim["/PartType0/InternalEnergies"][:]
# Compute the temperature
u *= unit_length_in_si ** 2 / unit_time_in_si ** 2
diff --git a/examples/SantaBarbara/SantaBarbara-256/rhoTPlot.py b/examples/SantaBarbara/SantaBarbara-256/rhoTPlot.py
index c290268eaa548e188bb652104ea9e726ea88a267..3bcf01d2a49bc1c53b243ffcff12359201d26d87 100644
--- a/examples/SantaBarbara/SantaBarbara-256/rhoTPlot.py
+++ b/examples/SantaBarbara/SantaBarbara-256/rhoTPlot.py
@@ -28,10 +28,10 @@ def get_data(filename):
data = SWIFTDataset(filename)
- data.gas.density.convert_to_units(mh / (cm ** 3))
- data.gas.temperature.convert_to_cgs()
+ data.gas.densities.convert_to_units(mh / (cm ** 3))
+ data.gas.temperatures.convert_to_cgs()
- return data.gas.density, data.gas.temperature
+ return data.gas.densities, data.gas.temperatures
def make_hist(filename, density_bounds, temperature_bounds, bins):
@@ -155,10 +155,8 @@ def make_movie(args, density_bounds, temperature_bounds, bins):
def format_metadata(metadata: SWIFTMetadata):
t = metadata.t * units.units["Unit time in cgs (U_t)"]
t.convert_to_units(Gyr)
-
- x = "$a$: {:2.2f}\n$z$: {:2.2f}\n$t$: {:2.2f}".format(
- metadata.a, metadata.z, t
- )
+
+ x = "$a$: {:2.2f}\n$z$: {:2.2f}\n$t$: {:2.2f}".format(metadata.a, metadata.z, t)
return x
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotRhoT.py b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotRhoT.py
index 4ba8ad66daca1d9614be8917a77407dd99209dea..4f02213ec2a66700d28ad5f8e57e00c30f3019d7 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotRhoT.py
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotRhoT.py
@@ -112,8 +112,8 @@ def T(u, H_frac=H_mass_fraction, T_trans=H_transition_temp):
return ret
-rho = sim["/PartType0/Density"][:]
-u = sim["/PartType0/InternalEnergy"][:]
+rho = sim["/PartType0/Densities"][:]
+u = sim["/PartType0/InternalEnergies"][:]
# Compute the temperature
u *= unit_length_in_si ** 2 / unit_time_in_si ** 2
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotTempEvolution.py b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotTempEvolution.py
index a3458ac1598e5657f3f597dfb10b36a7a641e68f..1e8cf9ea1082372d8e395c352f908c7ce693d99f 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotTempEvolution.py
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotTempEvolution.py
@@ -126,7 +126,7 @@ for i in range(n_snapshots):
z[i] = sim["/Cosmology"].attrs["Redshift"][0]
a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
- u = sim["/PartType0/InternalEnergy"][:]
+ u = sim["/PartType0/InternalEnergies"][:]
# Compute the temperature
u *= (unit_length_in_si**2 / unit_time_in_si**2)
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/plotTempEvolution.py b/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/plotTempEvolution.py
index aa6c5df5fe5ff5c7d0944a45bb11344f70c57844..d707f70450471f2d2fc589dbc382366280e0e7f3 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/plotTempEvolution.py
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/plotTempEvolution.py
@@ -123,7 +123,7 @@ for i in range(n_snapshots):
z[i] = sim["/Cosmology"].attrs["Redshift"][0]
a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
- u = sim["/PartType0/InternalEnergy"][:]
+ u = sim["/PartType0/InternalEnergies"][:]
# Compute the temperature
u *= (unit_length_in_si**2 / unit_time_in_si**2)
diff --git a/examples/SubgridTests/BlackHoleSwallowing/check_masses.py b/examples/SubgridTests/BlackHoleSwallowing/check_masses.py
index c7f1d7b2c1f90efa285c13c75f0e5243f36e49ea..a5b55ed00df0851f989858ddffd00ea34df88a27 100644
--- a/examples/SubgridTests/BlackHoleSwallowing/check_masses.py
+++ b/examples/SubgridTests/BlackHoleSwallowing/check_masses.py
@@ -66,5 +66,5 @@ for i in range(np.size(ids_removed)):
result = np.where(ids_gas == ids_removed)
print result
-#rho_gas = f["/PartType0/Density"][:]
+#rho_gas = f["/PartType0/Densities"][:]
#print np.mean(rho_gas), np.std(rho_gas)
diff --git a/examples/SubgridTests/SmoothedMetallicity/plotSolution.py b/examples/SubgridTests/SmoothedMetallicity/plotSolution.py
index e5bca3dfb7fe1e43c836733894c9e297cdd468ca..068fe5378e19c34ee8a68398f4e0ed096d0982e0 100644
--- a/examples/SubgridTests/SmoothedMetallicity/plotSolution.py
+++ b/examples/SubgridTests/SmoothedMetallicity/plotSolution.py
@@ -3,20 +3,20 @@
# This file is part of SWIFT.
# Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
# Matthieu Schaller (matthieu.schaller@durham.ac.uk)
-#
+#
# 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/>.
-#
+#
##############################################################################
# Computes the analytical solution of the 3D Smoothed Metallicity example.
@@ -25,12 +25,13 @@ import h5py
import sys
import numpy as np
import matplotlib
+
matplotlib.use("Agg")
import matplotlib.pyplot as plt
# Parameters
-low_metal = -6 # low metal abundance
-high_metal = -5 # High metal abundance
+low_metal = -6 # low metal abundance
+high_metal = -5 # High metal abundance
sigma_metal = 0.1 # relative standard deviation for Z
Nelem = 9
@@ -44,27 +45,27 @@ high_metal = [high_metal] * Nelem + np.linspace(0, 3, Nelem)
# Plot parameters
params = {
- 'axes.labelsize': 10,
- 'axes.titlesize': 10,
- 'font.size': 12,
- 'legend.fontsize': 12,
- 'xtick.labelsize': 10,
- 'ytick.labelsize': 10,
- 'text.usetex': True,
- 'figure.figsize': (9.90, 6.45),
- 'figure.subplot.left': 0.045,
- 'figure.subplot.right': 0.99,
- 'figure.subplot.bottom': 0.05,
- 'figure.subplot.top': 0.99,
- 'figure.subplot.wspace': 0.15,
- 'figure.subplot.hspace': 0.12,
- 'lines.markersize': 6,
- 'lines.linewidth': 3.,
- 'text.latex.unicode': True
+ "axes.labelsize": 10,
+ "axes.titlesize": 10,
+ "font.size": 12,
+ "legend.fontsize": 12,
+ "xtick.labelsize": 10,
+ "ytick.labelsize": 10,
+ "text.usetex": True,
+ "figure.figsize": (9.90, 6.45),
+ "figure.subplot.left": 0.045,
+ "figure.subplot.right": 0.99,
+ "figure.subplot.bottom": 0.05,
+ "figure.subplot.top": 0.99,
+ "figure.subplot.wspace": 0.15,
+ "figure.subplot.hspace": 0.12,
+ "lines.markersize": 6,
+ "lines.linewidth": 3.0,
+ "text.latex.unicode": True,
}
plt.rcParams.update(params)
-plt.rc('font', **{'family': 'sans-serif', 'sans-serif': ['Times']})
+plt.rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
snap = int(sys.argv[1])
@@ -83,18 +84,17 @@ git = sim["Code"].attrs["Git Revision"]
pos = sim["/PartType0/Coordinates"][:, :]
d = pos[:, 0] - boxSize / 2
-smooth_metal = sim["/PartType0/SmoothedElementAbundance"][:, :]
-metal = sim["/PartType0/ElementAbundance"][:, :]
-h = sim["/PartType0/SmoothingLength"][:]
+smooth_metal = sim["/PartType0/SmoothedElementMassFractions"][:, :]
+metal = sim["/PartType0/ElementMassFractions"][:, :]
+h = sim["/PartType0/SmoothingLengths"][:]
h = np.mean(h)
-if (Nelem != metal.shape[1]):
- print("Unexpected number of element, please check makeIC.py"
- " and plotSolution.py")
+if Nelem != metal.shape[1]:
+ print("Unexpected number of element, please check makeIC.py" " and plotSolution.py")
exit(1)
N = 1000
-d_a = np.linspace(-boxSize / 2., boxSize / 2., N)
+d_a = np.linspace(-boxSize / 2.0, boxSize / 2.0, N)
# Now, work our the solution....
@@ -142,14 +142,14 @@ def calc_a(d, high_metal, low_metal, std_dev, h):
m = (high_metal[i] - low_metal[i]) / (2.0 * h)
c = (high_metal[i] + low_metal[i]) / 2.0
# compute left linear part
- s = d < - boxSize / 2.0 + h
- a[s, i] = - m * (d[s] + boxSize / 2.0) + c
+ s = d < -boxSize / 2.0 + h
+ a[s, i] = -m * (d[s] + boxSize / 2.0) + c
# compute middle linear part
s = np.logical_and(d >= -h, d <= h)
a[s, i] = m * d[s] + c
# compute right linear part
s = d > boxSize / 2.0 - h
- a[s, i] = - m * (d[s] - boxSize / 2.0) + c
+ a[s, i] = -m * (d[s] - boxSize / 2.0) + c
sigma[:, :, 0] = a * (1 + std_dev)
sigma[:, :, 1] = a * (1 - std_dev)
@@ -165,7 +165,7 @@ plt.figure()
# Metallicity --------------------------------
plt.subplot(221)
for e in range(Nelem):
- plt.plot(metal[:, e], smooth_metal[:, e], '.', ms=0.5, alpha=0.2)
+ plt.plot(metal[:, e], smooth_metal[:, e], ".", ms=0.5, alpha=0.2)
xmin, xmax = metal.min(), metal.max()
ymin, ymax = smooth_metal.min(), smooth_metal.max()
@@ -178,27 +178,28 @@ plt.ylabel("${\\rm{Smoothed~Metallicity}}~Z_\\textrm{sm}$", labelpad=0)
# Metallicity --------------------------------
e = 0
plt.subplot(223)
-plt.plot(d, smooth_metal[:, e], '.', color='r', ms=0.5, alpha=0.2)
-plt.plot(d_a, sol[:, e], '--', color='b', alpha=0.8, lw=1.2)
-plt.fill_between(d_a, sig[:, e, 0], sig[:, e, 1], facecolor="b",
- interpolate=True, alpha=0.5)
+plt.plot(d, smooth_metal[:, e], ".", color="r", ms=0.5, alpha=0.2)
+plt.plot(d_a, sol[:, e], "--", color="b", alpha=0.8, lw=1.2)
+plt.fill_between(
+ d_a, sig[:, e, 0], sig[:, e, 1], facecolor="b", interpolate=True, alpha=0.5
+)
plt.xlabel("${\\rm{Distance}}~r$", labelpad=0)
plt.ylabel("${\\rm{Smoothed~Metallicity}}~Z_\\textrm{sm}$", labelpad=0)
plt.xlim(-0.5, 0.5)
-plt.ylim(low_metal[e]-1, high_metal[e]+1)
+plt.ylim(low_metal[e] - 1, high_metal[e] + 1)
# Information -------------------------------------
plt.subplot(222, frameon=False)
-plt.text(-0.49, 0.9, "Smoothed Metallicity in 3D at $t=%.2f$" % time,
- fontsize=10)
-plt.plot([-0.49, 0.1], [0.82, 0.82], 'k-', lw=1)
+plt.text(-0.49, 0.9, "Smoothed Metallicity in 3D at $t=%.2f$" % time, fontsize=10)
+plt.plot([-0.49, 0.1], [0.82, 0.82], "k-", lw=1)
plt.text(-0.49, 0.7, "$\\textsc{Swift}$ %s" % git, fontsize=10)
plt.text(-0.49, 0.6, scheme, fontsize=10)
plt.text(-0.49, 0.5, kernel, fontsize=10)
plt.text(-0.49, 0.4, chemistry + "'s Chemistry", fontsize=10)
-plt.text(-0.49, 0.3, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
- fontsize=10)
+plt.text(
+ -0.49, 0.3, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta), fontsize=10
+)
plt.xlim(-0.5, 0.5)
plt.ylim(0, 1)
plt.xticks([])
diff --git a/examples/SubgridTests/StellarEvolution/check_continuous_heating.py b/examples/SubgridTests/StellarEvolution/check_continuous_heating.py
index f3c1b5d7fd682d914f2dbc05259c2dab0baf1e32..b5940eba43c89b8af4a883cc8f7022e33293b869 100644
--- a/examples/SubgridTests/StellarEvolution/check_continuous_heating.py
+++ b/examples/SubgridTests/StellarEvolution/check_continuous_heating.py
@@ -93,7 +93,7 @@ for i in range(n_snapshots):
sim = h5py.File("stellar_evolution_%04d.hdf5"%i, "r")
print('reading snapshot '+str(i))
masses[:,i] = sim["/PartType0/Masses"]
- internal_energy[:,i] = sim["/PartType0/InternalEnergy"]
+ internal_energy[:,i] = sim["/PartType0/InternalEnergies"]
velocity_parts[:,:,i] = sim["/PartType0/Velocities"]
time[i] = sim["/Header"].attrs["Time"][0]
diff --git a/examples/SubgridTests/StellarEvolution/check_stellar_evolution.py b/examples/SubgridTests/StellarEvolution/check_stellar_evolution.py
index 02c1e9343de7b58cfddc8dee3bf0215a4b80ccf4..5680eb4d64f29ab32831d995e31ae9c27de82a71 100644
--- a/examples/SubgridTests/StellarEvolution/check_stellar_evolution.py
+++ b/examples/SubgridTests/StellarEvolution/check_stellar_evolution.py
@@ -1,4 +1,5 @@
import matplotlib
+
matplotlib.use("Agg")
from pylab import *
import h5py
@@ -7,32 +8,35 @@ import numpy as np
import glob
# Number of snapshots and elements
-newest_snap_name = max(glob.glob('stellar_evolution_*.hdf5'), key=os.path.getctime)
-n_snapshots = int(newest_snap_name.replace('stellar_evolution_','').replace('.hdf5','')) + 1
+newest_snap_name = max(glob.glob("stellar_evolution_*.hdf5"), key=os.path.getctime)
+n_snapshots = (
+ int(newest_snap_name.replace("stellar_evolution_", "").replace(".hdf5", "")) + 1
+)
n_elements = 9
# Plot parameters
-params = {'axes.labelsize': 10,
-'axes.titlesize': 10,
-'font.size': 9,
-'legend.fontsize': 9,
-'xtick.labelsize': 10,
-'ytick.labelsize': 10,
-'text.usetex': True,
- 'figure.figsize' : (3.15,3.15),
-'figure.subplot.left' : 0.3,
-'figure.subplot.right' : 0.99,
-'figure.subplot.bottom' : 0.18,
-'figure.subplot.top' : 0.92,
-'figure.subplot.wspace' : 0.21,
-'figure.subplot.hspace' : 0.19,
-'lines.markersize' : 6,
-'lines.linewidth' : 2.,
-'text.latex.unicode': True
+params = {
+ "axes.labelsize": 10,
+ "axes.titlesize": 10,
+ "font.size": 9,
+ "legend.fontsize": 9,
+ "xtick.labelsize": 10,
+ "ytick.labelsize": 10,
+ "text.usetex": True,
+ "figure.figsize": (3.15, 3.15),
+ "figure.subplot.left": 0.3,
+ "figure.subplot.right": 0.99,
+ "figure.subplot.bottom": 0.18,
+ "figure.subplot.top": 0.92,
+ "figure.subplot.wspace": 0.21,
+ "figure.subplot.hspace": 0.19,
+ "lines.markersize": 6,
+ "lines.linewidth": 2.0,
+ "text.latex.unicode": True,
}
rcParams.update(params)
-rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
# Read the simulation data
sim = h5py.File("stellar_evolution_0000.hdf5", "r")
@@ -43,7 +47,9 @@ neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"][0]
eta = sim["/HydroScheme"].attrs["Kernel eta"][0]
alpha = sim["/HydroScheme"].attrs["Alpha viscosity"][0]
H_mass_fraction = sim["/HydroScheme"].attrs["Hydrogen mass fraction"][0]
-H_transition_temp = sim["/HydroScheme"].attrs["Hydrogen ionization transition temperature"][0]
+H_transition_temp = sim["/HydroScheme"].attrs[
+ "Hydrogen ionization transition temperature"
+][0]
T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
git = sim["Code"].attrs["Git Revision"]
@@ -68,136 +74,245 @@ n_parts = sim["/Header"].attrs["NumPart_Total"][0]
n_sparts = sim["/Header"].attrs["NumPart_Total"][4]
# Declare arrays for data
-masses = zeros((n_parts,n_snapshots))
-star_masses = zeros((n_sparts,n_snapshots))
-mass_from_AGB = zeros((n_parts,n_snapshots))
-metal_mass_frac_from_AGB = zeros((n_parts,n_snapshots))
-mass_from_SNII = zeros((n_parts,n_snapshots))
-metal_mass_frac_from_SNII = zeros((n_parts,n_snapshots))
-mass_from_SNIa = zeros((n_parts,n_snapshots))
-metal_mass_frac_from_SNIa = zeros((n_parts,n_snapshots))
-iron_mass_frac_from_SNIa = zeros((n_parts,n_snapshots))
-metallicity = zeros((n_parts,n_snapshots))
-abundances = zeros((n_parts,n_elements,n_snapshots))
-internal_energy = zeros((n_parts,n_snapshots))
-coord_parts = zeros((n_parts,3))
-velocity_parts = zeros((n_parts,3,n_snapshots))
+masses = zeros((n_parts, n_snapshots))
+star_masses = zeros((n_sparts, n_snapshots))
+mass_from_AGB = zeros((n_parts, n_snapshots))
+metal_mass_frac_from_AGB = zeros((n_parts, n_snapshots))
+mass_from_SNII = zeros((n_parts, n_snapshots))
+metal_mass_frac_from_SNII = zeros((n_parts, n_snapshots))
+mass_from_SNIa = zeros((n_parts, n_snapshots))
+metal_mass_frac_from_SNIa = zeros((n_parts, n_snapshots))
+iron_mass_frac_from_SNIa = zeros((n_parts, n_snapshots))
+metallicity = zeros((n_parts, n_snapshots))
+abundances = zeros((n_parts, n_elements, n_snapshots))
+internal_energy = zeros((n_parts, n_snapshots))
+coord_parts = zeros((n_parts, 3))
+velocity_parts = zeros((n_parts, 3, n_snapshots))
coord_sparts = zeros(3)
time = zeros(n_snapshots)
# Read fields we are checking from snapshots
for i in range(n_snapshots):
- sim = h5py.File("stellar_evolution_%04d.hdf5"%i, "r")
- print('reading snapshot '+str(i))
- abundances[:,:,i] = sim["/PartType0/ElementAbundance"]
- metallicity[:,i] = sim["/PartType0/Metallicity"]
- masses[:,i] = sim["/PartType0/Masses"]
- star_masses[:,i] = sim["/PartType4/Masses"]
- mass_from_AGB[:,i] = sim["/PartType0/TotalMassFromAGB"]
- metal_mass_frac_from_AGB[:,i] = sim["/PartType0/MetalMassFracFromAGB"]
- mass_from_SNII[:,i] = sim["/PartType0/TotalMassFromSNII"]
- metal_mass_frac_from_SNII[:,i] = sim["/PartType0/MetalMassFracFromSNII"]
- mass_from_SNIa[:,i] = sim["/PartType0/TotalMassFromSNIa"]
- metal_mass_frac_from_SNIa[:,i] = sim["/PartType0/MetalMassFracFromSNIa"]
- iron_mass_frac_from_SNIa[:,i] = sim["/PartType0/IronMassFracFromSNIa"]
- internal_energy[:,i] = sim["/PartType0/InternalEnergy"]
- velocity_parts[:,:,i] = sim["/PartType0/Velocities"]
- time[i] = sim["/Header"].attrs["Time"][0]
+ sim = h5py.File("stellar_evolution_%04d.hdf5" % i, "r")
+ print("reading snapshot " + str(i))
+ abundances[:, :, i] = sim["/PartType0/ElementMassFractions"]
+ metallicity[:, i] = sim["/PartType0/Metallicities"]
+ masses[:, i] = sim["/PartType0/Masses"]
+ star_masses[:, i] = sim["/PartType4/Masses"]
+ mass_from_AGB[:, i] = sim["/PartType0/TotalMassFromAGB"]
+ metal_mass_frac_from_AGB[:, i] = sim["/PartType0/MetalMassFracFromAGB"]
+ mass_from_SNII[:, i] = sim["/PartType0/TotalMassFromSNII"]
+ metal_mass_frac_from_SNII[:, i] = sim["/PartType0/MetalMassFracFromSNII"]
+ mass_from_SNIa[:, i] = sim["/PartType0/TotalMassFromSNIa"]
+ metal_mass_frac_from_SNIa[:, i] = sim["/PartType0/MetalMassFracFromSNIa"]
+ iron_mass_frac_from_SNIa[:, i] = sim["/PartType0/IronMassFracFromSNIa"]
+ internal_energy[:, i] = sim["/PartType0/InternalEnergies"]
+ velocity_parts[:, :, i] = sim["/PartType0/Velocities"]
+ time[i] = sim["/Header"].attrs["Time"][0]
# Define ejecta factor
ejecta_factor = 1.0e-2
-ejecta_factor_metallicity = 1.0 - 2.0/n_elements
-ejecta_factor_abundances = 1.0/n_elements
+ejecta_factor_metallicity = 1.0 - 2.0 / n_elements
+ejecta_factor_abundances = 1.0 / n_elements
ejected_mass = star_initial_mass
-energy_per_SNe = 1.0e51/unit_energy_in_cgs
+energy_per_SNe = 1.0e51 / unit_energy_in_cgs
# Check that the total amount of enrichment is as expected.
# Define tolerance
eps = 0.01
# Total mass
-total_part_mass = np.sum(masses,axis = 0)
-if abs((total_part_mass[n_snapshots-1] - total_part_mass[0])/total_part_mass[0] - ejected_mass/total_part_mass[0])*total_part_mass[0]/ejected_mass < eps:
- print("total mass released consistent with expectation")
+total_part_mass = np.sum(masses, axis=0)
+if (
+ abs(
+ (total_part_mass[n_snapshots - 1] - total_part_mass[0]) / total_part_mass[0]
+ - ejected_mass / total_part_mass[0]
+ )
+ * total_part_mass[0]
+ / ejected_mass
+ < eps
+):
+ print("total mass released consistent with expectation")
else:
- print("mass increase "+str(total_part_mass[n_snapshots-1]/total_part_mass[0])+" expected "+ str(1.0+ejected_mass/total_part_mass[0]))
+ print(
+ "mass increase "
+ + str(total_part_mass[n_snapshots - 1] / total_part_mass[0])
+ + " expected "
+ + str(1.0 + ejected_mass / total_part_mass[0])
+ )
# Check that mass is conserved (i.e. total star mass decreases by same amount as total gas mass increases)
-total_spart_mass = np.sum(star_masses,axis = 0)
-if abs((total_part_mass[n_snapshots-1] + total_spart_mass[n_snapshots-1]) / (total_part_mass[0] + total_spart_mass[0]) - 1.0) < eps**3:
- print("total mass conserved")
+total_spart_mass = np.sum(star_masses, axis=0)
+if (
+ abs(
+ (total_part_mass[n_snapshots - 1] + total_spart_mass[n_snapshots - 1])
+ / (total_part_mass[0] + total_spart_mass[0])
+ - 1.0
+ )
+ < eps ** 3
+):
+ print("total mass conserved")
else:
- print("initial part, spart mass " + str(total_part_mass[0]) + " " + str(total_spart_mass[0]) + " final mass " + str(total_part_mass[n_snapshots-1]) + " " + str(total_spart_mass[n_snapshots-1]))
+ print(
+ "initial part, spart mass "
+ + str(total_part_mass[0])
+ + " "
+ + str(total_spart_mass[0])
+ + " final mass "
+ + str(total_part_mass[n_snapshots - 1])
+ + " "
+ + str(total_spart_mass[n_snapshots - 1])
+ )
# Total metal mass from AGB
-total_metal_mass_AGB = np.sum(np.multiply(metal_mass_frac_from_AGB,masses),axis = 0)
-expected_metal_mass_AGB = ejecta_factor*ejected_mass
-if abs(total_metal_mass_AGB[n_snapshots-1] - expected_metal_mass_AGB)/expected_metal_mass_AGB < eps:
- print("total AGB metal mass released consistent with expectation")
+total_metal_mass_AGB = np.sum(np.multiply(metal_mass_frac_from_AGB, masses), axis=0)
+expected_metal_mass_AGB = ejecta_factor * ejected_mass
+if (
+ abs(total_metal_mass_AGB[n_snapshots - 1] - expected_metal_mass_AGB)
+ / expected_metal_mass_AGB
+ < eps
+):
+ print("total AGB metal mass released consistent with expectation")
else:
- print("total AGB metal mass "+str(total_metal_mass_AGB[n_snapshots-1])+" expected "+ str(expected_metal_mass_AGB))
+ print(
+ "total AGB metal mass "
+ + str(total_metal_mass_AGB[n_snapshots - 1])
+ + " expected "
+ + str(expected_metal_mass_AGB)
+ )
# Total mass from AGB
-total_AGB_mass = np.sum(mass_from_AGB,axis = 0)
-expected_AGB_mass = ejecta_factor*ejected_mass
-if abs(total_AGB_mass[n_snapshots-1] - expected_AGB_mass)/expected_AGB_mass < eps:
- print("total AGB mass released consistent with expectation")
+total_AGB_mass = np.sum(mass_from_AGB, axis=0)
+expected_AGB_mass = ejecta_factor * ejected_mass
+if abs(total_AGB_mass[n_snapshots - 1] - expected_AGB_mass) / expected_AGB_mass < eps:
+ print("total AGB mass released consistent with expectation")
else:
- print("total AGB mass "+str(total_AGB_mass[n_snapshots-1])+" expected "+ str(expected_AGB_mass))
+ print(
+ "total AGB mass "
+ + str(total_AGB_mass[n_snapshots - 1])
+ + " expected "
+ + str(expected_AGB_mass)
+ )
# Total metal mass from SNII
-total_metal_mass_SNII = np.sum(np.multiply(metal_mass_frac_from_SNII,masses),axis = 0)
-expected_metal_mass_SNII = ejecta_factor*ejected_mass
-if abs(total_metal_mass_SNII[n_snapshots-1] - expected_metal_mass_SNII)/expected_metal_mass_SNII < eps:
- print("total SNII metal mass released consistent with expectation")
+total_metal_mass_SNII = np.sum(np.multiply(metal_mass_frac_from_SNII, masses), axis=0)
+expected_metal_mass_SNII = ejecta_factor * ejected_mass
+if (
+ abs(total_metal_mass_SNII[n_snapshots - 1] - expected_metal_mass_SNII)
+ / expected_metal_mass_SNII
+ < eps
+):
+ print("total SNII metal mass released consistent with expectation")
else:
- print("total SNII metal mass "+str(total_metal_mass_SNII[n_snapshots-1])+" expected "+ str(expected_metal_mass_SNII))
+ print(
+ "total SNII metal mass "
+ + str(total_metal_mass_SNII[n_snapshots - 1])
+ + " expected "
+ + str(expected_metal_mass_SNII)
+ )
# Total mass from SNII
-total_SNII_mass = np.sum(mass_from_SNII,axis = 0)
-expected_SNII_mass = ejecta_factor*ejected_mass
-if abs(total_SNII_mass[n_snapshots-1] - expected_SNII_mass)/expected_SNII_mass < eps:
- print("total SNII mass released consistent with expectation")
+total_SNII_mass = np.sum(mass_from_SNII, axis=0)
+expected_SNII_mass = ejecta_factor * ejected_mass
+if (
+ abs(total_SNII_mass[n_snapshots - 1] - expected_SNII_mass) / expected_SNII_mass
+ < eps
+):
+ print("total SNII mass released consistent with expectation")
else:
- print("total SNII mass "+str(total_SNII_mass[n_snapshots-1])+" expected "+ str(expected_SNII_mass))
+ print(
+ "total SNII mass "
+ + str(total_SNII_mass[n_snapshots - 1])
+ + " expected "
+ + str(expected_SNII_mass)
+ )
# Total metal mass from SNIa
-total_metal_mass_SNIa = np.sum(np.multiply(metal_mass_frac_from_SNIa,masses),axis = 0)
-expected_metal_mass_SNIa = ejecta_factor*ejected_mass
-if abs(total_metal_mass_SNIa[n_snapshots-1] - expected_metal_mass_SNIa)/expected_metal_mass_SNIa < eps:
- print("total SNIa metal mass released consistent with expectation")
+total_metal_mass_SNIa = np.sum(np.multiply(metal_mass_frac_from_SNIa, masses), axis=0)
+expected_metal_mass_SNIa = ejecta_factor * ejected_mass
+if (
+ abs(total_metal_mass_SNIa[n_snapshots - 1] - expected_metal_mass_SNIa)
+ / expected_metal_mass_SNIa
+ < eps
+):
+ print("total SNIa metal mass released consistent with expectation")
else:
- print("total SNIa metal mass "+str(total_metal_mass_SNIa[n_snapshots-1])+" expected "+ str(expected_metal_mass_SNIa))
+ print(
+ "total SNIa metal mass "
+ + str(total_metal_mass_SNIa[n_snapshots - 1])
+ + " expected "
+ + str(expected_metal_mass_SNIa)
+ )
# Total iron mass from SNIa
-total_iron_mass_SNIa = np.sum(np.multiply(iron_mass_frac_from_SNIa,masses),axis = 0)
-expected_iron_mass_SNIa = ejecta_factor*ejected_mass
-if abs(total_iron_mass_SNIa[n_snapshots-1] - expected_iron_mass_SNIa)/expected_iron_mass_SNIa < eps:
- print("total SNIa iron mass released consistent with expectation")
+total_iron_mass_SNIa = np.sum(np.multiply(iron_mass_frac_from_SNIa, masses), axis=0)
+expected_iron_mass_SNIa = ejecta_factor * ejected_mass
+if (
+ abs(total_iron_mass_SNIa[n_snapshots - 1] - expected_iron_mass_SNIa)
+ / expected_iron_mass_SNIa
+ < eps
+):
+ print("total SNIa iron mass released consistent with expectation")
else:
- print("total SNIa iron mass "+str(total_iron_mass_SNIa[n_snapshots-1])+" expected "+ str(expected_iron_mass_SNIa))
+ print(
+ "total SNIa iron mass "
+ + str(total_iron_mass_SNIa[n_snapshots - 1])
+ + " expected "
+ + str(expected_iron_mass_SNIa)
+ )
# Total mass from SNIa
-total_SNIa_mass = np.sum(mass_from_SNIa,axis = 0)
-expected_SNIa_mass = ejecta_factor*ejected_mass
-if abs(total_SNIa_mass[n_snapshots-1] - expected_SNIa_mass)/expected_SNIa_mass < eps:
- print("total SNIa mass released consistent with expectation")
+total_SNIa_mass = np.sum(mass_from_SNIa, axis=0)
+expected_SNIa_mass = ejecta_factor * ejected_mass
+if (
+ abs(total_SNIa_mass[n_snapshots - 1] - expected_SNIa_mass) / expected_SNIa_mass
+ < eps
+):
+ print("total SNIa mass released consistent with expectation")
else:
- print("total SNIa mass "+str(total_SNIa_mass[n_snapshots-1])+" expected "+ str(expected_SNIa_mass))
+ print(
+ "total SNIa mass "
+ + str(total_SNIa_mass[n_snapshots - 1])
+ + " expected "
+ + str(expected_SNIa_mass)
+ )
# Total metal mass
-total_metal_mass = np.sum(np.multiply(metallicity,masses),axis = 0)
-expected_metal_mass = ejecta_factor_metallicity*ejected_mass
-if abs(total_metal_mass[n_snapshots-1] - expected_metal_mass)/expected_metal_mass < eps:
- print("total metal mass released consistent with expectation")
+total_metal_mass = np.sum(np.multiply(metallicity, masses), axis=0)
+expected_metal_mass = ejecta_factor_metallicity * ejected_mass
+if (
+ abs(total_metal_mass[n_snapshots - 1] - expected_metal_mass) / expected_metal_mass
+ < eps
+):
+ print("total metal mass released consistent with expectation")
else:
- print("total metal mass "+str(total_metal_mass[n_snapshots-1])+" expected "+ str(expected_metal_mass))
+ print(
+ "total metal mass "
+ + str(total_metal_mass[n_snapshots - 1])
+ + " expected "
+ + str(expected_metal_mass)
+ )
# Total mass for each element
-expected_element_mass = ejecta_factor_abundances*ejected_mass
+expected_element_mass = ejecta_factor_abundances * ejected_mass
for i in range(n_elements):
- total_element_mass = np.sum(np.multiply(abundances[:,i,:],masses),axis = 0)
- if abs(total_element_mass[n_snapshots-1] - expected_element_mass)/expected_element_mass < eps:
- print("total element mass released consistent with expectation for element "+str(i))
- else:
- print("total element mass "+str(total_element_mass[n_snapshots-1])+" expected "+ str(expected_element_mass) + " for element "+ str(i))
+ total_element_mass = np.sum(np.multiply(abundances[:, i, :], masses), axis=0)
+ if (
+ abs(total_element_mass[n_snapshots - 1] - expected_element_mass)
+ / expected_element_mass
+ < eps
+ ):
+ print(
+ "total element mass released consistent with expectation for element "
+ + str(i)
+ )
+ else:
+ print(
+ "total element mass "
+ + str(total_element_mass[n_snapshots - 1])
+ + " expected "
+ + str(expected_element_mass)
+ + " for element "
+ + str(i)
+ )
+
diff --git a/examples/SubgridTests/StellarEvolution/check_stochastic_heating.py b/examples/SubgridTests/StellarEvolution/check_stochastic_heating.py
index da837540041a9295a33b55e16b5e996394576cd7..1cacc13653d821da3abd2a09566be347608c64f7 100644
--- a/examples/SubgridTests/StellarEvolution/check_stochastic_heating.py
+++ b/examples/SubgridTests/StellarEvolution/check_stochastic_heating.py
@@ -93,7 +93,7 @@ for i in range(n_snapshots):
sim = h5py.File("stellar_evolution_%04d.hdf5"%i, "r")
print('reading snapshot '+str(i))
masses[:,i] = sim["/PartType0/Masses"]
- internal_energy[:,i] = sim["/PartType0/InternalEnergy"]
+ internal_energy[:,i] = sim["/PartType0/InternalEnergies"]
velocity_parts[:,:,i] = sim["/PartType0/Velocities"]
time[i] = sim["/Header"].attrs["Time"][0]
diff --git a/examples/SubgridTests/StellarEvolution/plot_box_evolution.py b/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
index a46db721e153ade2d386a5790e26a290cead4f90..bcfa85a1afac021e280a797641dd677bccd291d3 100644
--- a/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
+++ b/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
@@ -112,16 +112,16 @@ for i in range(n_snapshots):
star_masses = sim["/PartType4/Masses"][:]
swift_box_star_mass[i] = np.sum(star_masses)
- metallicities = sim["/PartType0/Metallicity"][:]
+ metallicities = sim["/PartType0/Metallicities"][:]
swift_box_gas_metal_mass[i] = np.sum(metallicities * masses)
- element_abundances = sim["/PartType0/ElementAbundance"][:][:]
+ element_abundances = sim["/PartType0/ElementMassFractions"][:][:]
for j in range(n_elements):
swift_element_mass[i,j] = np.sum(element_abundances[:,j] * masses)
v = sim["/PartType0/Velocities"][:,:]
v2 = v[:,0]**2 + v[:,1]**2 + v[:,2]**2
- u = sim["/PartType0/InternalEnergy"][:]
+ u = sim["/PartType0/InternalEnergies"][:]
swift_internal_energy[i] = np.sum(masses * u)
swift_kinetic_energy[i] = np.sum(0.5 * masses * v2)
swift_total_energy[i] = swift_kinetic_energy[i] + swift_internal_energy[i]
diff --git a/examples/SubgridTests/StellarEvolution/plot_particle_evolution.py b/examples/SubgridTests/StellarEvolution/plot_particle_evolution.py
index 8b935e537b14a9d1d9cc4eec7c5cd0794c6fc489..be1588f9b707d448b2905611defd9e760c5f91de 100644
--- a/examples/SubgridTests/StellarEvolution/plot_particle_evolution.py
+++ b/examples/SubgridTests/StellarEvolution/plot_particle_evolution.py
@@ -1,29 +1,30 @@
###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
- # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- #
- # 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/>.
- #
- ##############################################################################
-
-# Assuming output snapshots contain evolution of box of gas with star at its
-# centre, this script will plot the evolution of the radial velocities, internal
-# energies, mass and metallicities of the nearest n particles to the star over
-# the duration of the simulation.
+# This file is part of SWIFT.
+# Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+# Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#
+# 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/>.
+#
+##############################################################################
+
+# Assuming output snapshots contain evolution of box of gas with star at its
+# centre, this script will plot the evolution of the radial velocities, internal
+# energies, mass and metallicities of the nearest n particles to the star over
+# the duration of the simulation.
import matplotlib
+
matplotlib.use("Agg")
from pylab import *
from scipy import stats
@@ -33,38 +34,42 @@ import glob
import os.path
# Function to find index in array a for each element in array b
-def find_indices(a,b):
- result = np.zeros(len(b))
- for i in range(len(b)):
- result[i] = ((np.where(a == b[i]))[0])[0]
- return result
+def find_indices(a, b):
+ result = np.zeros(len(b))
+ for i in range(len(b)):
+ result[i] = ((np.where(a == b[i]))[0])[0]
+ return result
+
# Plot parameters
-params = {'axes.labelsize': 10,
-'axes.titlesize': 10,
-'font.size': 12,
-'legend.fontsize': 12,
-'xtick.labelsize': 10,
-'ytick.labelsize': 10,
-'text.usetex': True,
- 'figure.figsize' : (9.90,6.45),
-'figure.subplot.left' : 0.1,
-'figure.subplot.right' : 0.99,
-'figure.subplot.bottom' : 0.1,
-'figure.subplot.top' : 0.95,
-'figure.subplot.wspace' : 0.2,
-'figure.subplot.hspace' : 0.2,
-'lines.markersize' : 6,
-'lines.linewidth' : 3.,
-'text.latex.unicode': True
+params = {
+ "axes.labelsize": 10,
+ "axes.titlesize": 10,
+ "font.size": 12,
+ "legend.fontsize": 12,
+ "xtick.labelsize": 10,
+ "ytick.labelsize": 10,
+ "text.usetex": True,
+ "figure.figsize": (9.90, 6.45),
+ "figure.subplot.left": 0.1,
+ "figure.subplot.right": 0.99,
+ "figure.subplot.bottom": 0.1,
+ "figure.subplot.top": 0.95,
+ "figure.subplot.wspace": 0.2,
+ "figure.subplot.hspace": 0.2,
+ "lines.markersize": 6,
+ "lines.linewidth": 3.0,
+ "text.latex.unicode": True,
}
rcParams.update(params)
-rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
# Number of snapshots and elements
-newest_snap_name = max(glob.glob('stellar_evolution_*.hdf5'), key=os.path.getctime)
-n_snapshots = int(newest_snap_name.replace('stellar_evolution_','').replace('.hdf5','')) + 1
+newest_snap_name = max(glob.glob("stellar_evolution_*.hdf5"), key=os.path.getctime)
+n_snapshots = (
+ int(newest_snap_name.replace("stellar_evolution_", "").replace(".hdf5", "")) + 1
+)
n_particles_to_plot = 500
# Read the simulation data
@@ -87,25 +92,25 @@ unit_energy_in_cgs = unit_mass_in_cgs * unit_vel_in_cgs * unit_vel_in_cgs
unit_length_in_si = 0.01 * unit_length_in_cgs
unit_mass_in_si = 0.001 * unit_mass_in_cgs
unit_time_in_si = unit_time_in_cgs
-unit_density_in_cgs = unit_mass_in_cgs*unit_length_in_cgs**-3
-unit_pressure_in_cgs = unit_mass_in_cgs/unit_length_in_cgs*unit_time_in_cgs**-2
-unit_int_energy_in_cgs = unit_energy_in_cgs/unit_mass_in_cgs
-unit_entropy_in_cgs = unit_energy_in_cgs/unit_temp_in_cgs
+unit_density_in_cgs = unit_mass_in_cgs * unit_length_in_cgs ** -3
+unit_pressure_in_cgs = unit_mass_in_cgs / unit_length_in_cgs * unit_time_in_cgs ** -2
+unit_int_energy_in_cgs = unit_energy_in_cgs / unit_mass_in_cgs
+unit_entropy_in_cgs = unit_energy_in_cgs / unit_temp_in_cgs
Myr_in_cgs = 3.154e13
Msun_in_cgs = 1.989e33
# Read data of zeroth snapshot
-pos = sim["/PartType0/Coordinates"][:,:]
-x = pos[:,0] - boxSize / 2
-y = pos[:,1] - boxSize / 2
-z = pos[:,2] - boxSize / 2
-vel = sim["/PartType0/Velocities"][:,:]
-r = sqrt(x**2 + y**2 + z**2)
-v_r = (x * vel[:,0] + y * vel[:,1] + z * vel[:,2]) / r
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+pos = sim["/PartType0/Coordinates"][:, :]
+x = pos[:, 0] - boxSize / 2
+y = pos[:, 1] - boxSize / 2
+z = pos[:, 2] - boxSize / 2
+vel = sim["/PartType0/Velocities"][:, :]
+r = sqrt(x ** 2 + y ** 2 + z ** 2)
+v_r = (x * vel[:, 0] + y * vel[:, 1] + z * vel[:, 2]) / r
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
mass = sim["/PartType0/Masses"][:]
IDs = sim["/PartType0/ParticleIDs"][:]
@@ -123,34 +128,34 @@ t = zeros(n_snapshots)
# Read data from rest of snapshots
for i in range(n_snapshots):
- print("reading snapshot "+str(i))
- # Read the simulation data
- sim = h5py.File("stellar_evolution_%04d.hdf5"%i, "r")
- t[i] = sim["/Header"].attrs["Time"][0]
-
- pos = sim["/PartType0/Coordinates"][:,:]
- x = pos[:,0] - boxSize / 2
- y = pos[:,1] - boxSize / 2
- z = pos[:,2] - boxSize / 2
- vel = sim["/PartType0/Velocities"][:,:]
- r = sqrt(x**2 + y**2 + z**2)
- v_r = (x * vel[:,0] + y * vel[:,1] + z * vel[:,2]) / r
- u = sim["/PartType0/InternalEnergy"][:]
- S = sim["/PartType0/Entropy"][:]
- P = sim["/PartType0/Pressure"][:]
- rho = sim["/PartType0/Density"][:]
- mass = sim["/PartType0/Masses"][:]
- metallicity = sim["/PartType0/Metallicity"][:]
- internal_energy = sim["/PartType0/InternalEnergy"][:]
- IDs = sim["/PartType0/ParticleIDs"][:]
-
- # Find which particles we want to plot and store their data
- indices = (find_indices(IDs,part_IDs_to_plot)).astype(int)
- masses_to_plot[:,i] = mass[indices[:]]
- v_r_to_plot[:,i] = v_r[indices[:]]
- metallicities_to_plot[:,i] = metallicity[indices[:]]
- internal_energies_to_plot[:,i] = internal_energy[indices[:]]
-
+ print("reading snapshot " + str(i))
+ # Read the simulation data
+ sim = h5py.File("stellar_evolution_%04d.hdf5" % i, "r")
+ t[i] = sim["/Header"].attrs["Time"][0]
+
+ pos = sim["/PartType0/Coordinates"][:, :]
+ x = pos[:, 0] - boxSize / 2
+ y = pos[:, 1] - boxSize / 2
+ z = pos[:, 2] - boxSize / 2
+ vel = sim["/PartType0/Velocities"][:, :]
+ r = sqrt(x ** 2 + y ** 2 + z ** 2)
+ v_r = (x * vel[:, 0] + y * vel[:, 1] + z * vel[:, 2]) / r
+ u = sim["/PartType0/InternalEnergies"][:]
+ S = sim["/PartType0/Entropies"][:]
+ P = sim["/PartType0/Pressures"][:]
+ rho = sim["/PartType0/Densities"][:]
+ mass = sim["/PartType0/Masses"][:]
+ metallicity = sim["/PartType0/Metallicities"][:]
+ internal_energy = sim["/PartType0/InternalEnergies"][:]
+ IDs = sim["/PartType0/ParticleIDs"][:]
+
+ # Find which particles we want to plot and store their data
+ indices = (find_indices(IDs, part_IDs_to_plot)).astype(int)
+ masses_to_plot[:, i] = mass[indices[:]]
+ v_r_to_plot[:, i] = v_r[indices[:]]
+ metallicities_to_plot[:, i] = metallicity[indices[:]]
+ internal_energies_to_plot[:, i] = internal_energy[indices[:]]
+
# Plot the interesting quantities
figure()
@@ -158,33 +163,61 @@ figure()
# Radial velocity --------------------------------
subplot(221)
for j in range(n_particles_to_plot):
- plot(t * unit_time_in_cgs / Myr_in_cgs, v_r_to_plot[j,:] * unit_vel_in_cgs, linewidth=0.5, color='k', ms=0.5, alpha=0.1)
+ plot(
+ t * unit_time_in_cgs / Myr_in_cgs,
+ v_r_to_plot[j, :] * unit_vel_in_cgs,
+ linewidth=0.5,
+ color="k",
+ ms=0.5,
+ alpha=0.1,
+ )
xlabel("Time (Myr)", labelpad=0)
ylabel("Radial velocity $(\\rm{cm} \cdot \\rm{s}^{-1})$", labelpad=0)
-ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
# Internal energy --------------------------------
subplot(222)
for j in range(n_particles_to_plot):
- plot(t * unit_time_in_cgs / Myr_in_cgs, internal_energies_to_plot[j,:] * unit_energy_in_cgs / unit_mass_in_cgs, linewidth=0.5, color='k', ms=0.5, alpha=0.1)
+ plot(
+ t * unit_time_in_cgs / Myr_in_cgs,
+ internal_energies_to_plot[j, :] * unit_energy_in_cgs / unit_mass_in_cgs,
+ linewidth=0.5,
+ color="k",
+ ms=0.5,
+ alpha=0.1,
+ )
xlabel("Time (Myr)", labelpad=0)
ylabel("Internal energy $(\\rm{erg} \cdot \\rm{g}^{-1})$", labelpad=2)
-ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
# Masses --------------------------------
subplot(223)
for j in range(n_particles_to_plot):
- plot(t * unit_time_in_cgs / Myr_in_cgs, masses_to_plot[j,:] * unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', ms=0.5, alpha=0.1)
+ plot(
+ t * unit_time_in_cgs / Myr_in_cgs,
+ masses_to_plot[j, :] * unit_mass_in_cgs / Msun_in_cgs,
+ linewidth=0.5,
+ color="k",
+ ms=0.5,
+ alpha=0.1,
+ )
xlabel("Time (Myr)", labelpad=0)
ylabel("Mass (Msun)", labelpad=2)
-ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
# Metallicities --------------------------------
subplot(224)
for j in range(n_particles_to_plot):
- plot(t * unit_time_in_cgs / Myr_in_cgs, metallicities_to_plot[j,:] , linewidth=0.5, color='k', ms=0.5, alpha=0.1)
+ plot(
+ t * unit_time_in_cgs / Myr_in_cgs,
+ metallicities_to_plot[j, :],
+ linewidth=0.5,
+ color="k",
+ ms=0.5,
+ alpha=0.1,
+ )
xlabel("Time (Myr)", labelpad=0)
ylabel("Metallicity", labelpad=2)
-ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
savefig("particle_evolution.png", dpi=200)