diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/1Dprofiles.py b/examples/MHDTests/CoolingHaloWithSpin_MHD/1Dprofiles.py
index 7cdd0981a0f9fb6e4982e25d62650922d7393e80..ee568c98ee39df18e615f2c750589a68f2f81294 100644
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/1Dprofiles.py
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/1Dprofiles.py
@@ -11,14 +11,14 @@ from swiftsimio.visualisation.slice import slice_gas
PARSEC_IN_CGS = 3.0856776e18
KM_PER_SEC_IN_CGS = 1.0e5
CONST_G_CGS = 6.672e-8
-MSOL_IN_CGS = 1.9891e33 # Solar mass
-kb_cgs = 1.38e-16 # boltzmann constant
-m_H_cgs = 1.68e-24 # atomic hydrogen mass
-GYR_IN_CGS = 3.1536e16 # gigayear
+MSOL_IN_CGS = 1.9891e33 # Solar mass
+kb_cgs = 1.38e-16 # boltzmann constant
+m_H_cgs = 1.68e-24 # atomic hydrogen mass
+GYR_IN_CGS = 3.1536e16 # gigayear
# First set unit velocity and then the circular velocity parameter for the isothermal potential
const_unit_velocity_in_cgs = 1.0e5 # kms^-1
-h = 0.681 #0.67777 # hubble parameter
+h = 0.681 # 0.67777 # hubble parameter
H_0_cgs = 100.0 * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
# From this we can find the virial radius, the radius within which the average density of the halo is
@@ -29,12 +29,12 @@ f_b = 0.17
c_200 = 7.2
spin_lambda = 0.05
nH_max_cgs = 1e0
-M_200_cgs = 1e12 * MSOL_IN_CGS
-rhoc_cgs = 3*H_0_cgs**2/(8*np.pi*CONST_G_CGS)
-r_200_cgs = (3*M_200_cgs/(4*np.pi*rhoc_cgs*200))**(1/3)
-v_200_cgs = np.sqrt(CONST_G_CGS*M_200_cgs/r_200_cgs)
-v_200 = v_200_cgs / const_unit_velocity_in_cgs
-T_200_cgs = m_H_cgs*v_200_cgs**2/(2*kb_cgs)
+M_200_cgs = 1e12 * MSOL_IN_CGS
+rhoc_cgs = 3 * H_0_cgs ** 2 / (8 * np.pi * CONST_G_CGS)
+r_200_cgs = (3 * M_200_cgs / (4 * np.pi * rhoc_cgs * 200)) ** (1 / 3)
+v_200_cgs = np.sqrt(CONST_G_CGS * M_200_cgs / r_200_cgs)
+v_200 = v_200_cgs / const_unit_velocity_in_cgs
+T_200_cgs = m_H_cgs * v_200_cgs ** 2 / (2 * kb_cgs)
const_unit_mass_in_cgs = M_200_cgs
const_unit_length_in_cgs = r_200_cgs
@@ -48,15 +48,6 @@ const_G = (
/ (const_unit_length_in_cgs * const_unit_length_in_cgs * const_unit_length_in_cgs)
)
-# load reference
-with_reference = False
-
-if with_reference:
- import pandas as pd
- rho_vs_x_data = pd.read_csv('./1Dreference/rho_vs_x.csv',header=None)
- Bx_vs_x_data = pd.read_csv('./1Dreference/Bx_vs_x.csv',header=None)
- By_vs_x_data = pd.read_csv('./1Dreference/By_vs_x.csv',header=None)
- Berr_vs_x_data = pd.read_csv('./1Dreference/Berr_vs_x.csv',header=None)
# Parse command line arguments
argparser = argparse.ArgumentParser()
@@ -89,8 +80,10 @@ rho = data.gas.densities
m = data.gas.masses
coords = data.gas.coordinates
velocities = data.gas.velocities
-l = np.vstack([m,m,m]).T * np.cross((coords-data.metadata.boxsize/2),velocities) # specific angular momentum
-rho.convert_to_units(unyt.g*unyt.cm**(-3))
+l = np.vstack([m, m, m]).T * np.cross(
+ (coords - data.metadata.boxsize / 2), velocities
+) # specific angular momentum
+rho.convert_to_units(unyt.g * unyt.cm ** (-3))
P = data.gas.pressures
T = data.gas.temperatures
@@ -151,181 +144,231 @@ T_map = mass_weighted_T_map / mass_map
map_pixel_length = len(mass_map)
-n_H = density_map/m_H_cgs
+n_H = density_map / m_H_cgs
x = np.linspace(0.0, 1.0, map_pixel_length)
-slice_ind = int(np.round(y0 * map_pixel_length,0))
+slice_ind = int(np.round(y0 * map_pixel_length, 0))
# Plot maps
plt.rcParams.update({"font.size": 16})
nx = 1
-ny = 5
-fig, axs = plt.subplots(ny, nx, figsize=((10*nx, 5*ny)))
+ny = 4
+fig, axs = plt.subplots(ny, nx, figsize=((10 * nx, 5 * ny)))
fig.subplots_adjust(hspace=0.1)
-Lbox_kpc = data.metadata.boxsize.to(PARSEC_IN_CGS*1e3*unyt.cm).value[0]
-#locs = [map_pixel_length / 4, map_pixel_length / 2, 3 * map_pixel_length / 4]
-#labels = [-Lbox_kpc / 4, 0, Lbox_kpc / 4]
-x = np.linspace(-Lbox_kpc/2, Lbox_kpc/2, map_pixel_length)
+Lbox_kpc = data.metadata.boxsize.to(PARSEC_IN_CGS * 1e3 * unyt.cm).value[0]
+# locs = [map_pixel_length / 4, map_pixel_length / 2, 3 * map_pixel_length / 4]
+# labels = [-Lbox_kpc / 4, 0, Lbox_kpc / 4]
+x = np.linspace(-Lbox_kpc / 2, Lbox_kpc / 2, map_pixel_length)
slice_ind = int(np.floor(y0 * map_pixel_length))
# plot density
-axs[0].plot(x, n_H[:, slice_ind], "k-",color='black')
-axs[0].set_yscale('log')
+axs[0].plot(x, n_H[:, slice_ind], "k-", color="black")
+axs[0].set_yscale("log")
axs[0].set_yticks(np.logspace(-7, 2, 10))
axs[0].set_ylim(1e-7, 1e2)
# plot temperature
-P = pressure_map[:, slice_ind].to(MSOL_IN_CGS*unyt.g/(PARSEC_IN_CGS*unyt.cm * (GYR_IN_CGS*unyt.s)**2))
-axs[1].plot(x, P, "k-",color='black')
-axs[1].set_yscale('log')
+P = pressure_map[:, slice_ind].to(
+ MSOL_IN_CGS * unyt.g / (PARSEC_IN_CGS * unyt.cm * (GYR_IN_CGS * unyt.s) ** 2)
+)
+axs[1].plot(x, P, "k-", color="black")
+axs[1].set_yscale("log")
# plot mass
-coords = coords.to(1e3*PARSEC_IN_CGS*unyt.cm).value-Lbox_kpc/2
-r_values = np.logspace(np.log10(1e-6),np.log10(Lbox_kpc/2),100)
+coords = coords.to(1e3 * PARSEC_IN_CGS * unyt.cm).value - Lbox_kpc / 2
+r_values = np.logspace(np.log10(1e-6), np.log10(Lbox_kpc / 2), 100)
rp_kpc = np.linalg.norm(coords, axis=1)
-m=data.gas.masses.to(MSOL_IN_CGS*unyt.g).value
-M_x =np.array([ np.sum(m[rp_kpc<=r_values[i]]) for i in range(len(r_values))])
-axs[2].scatter(r_values,M_x,color='black', marker='+',s=100)
+m = data.gas.masses.to(MSOL_IN_CGS * unyt.g).value
+M_x = np.array([np.sum(m[rp_kpc <= r_values[i]]) for i in range(len(r_values))])
+axs[2].scatter(r_values, M_x, color="black", marker="+", s=100)
# plot specific angular momentum
-l_units_cgs = const_unit_mass_in_cgs * const_unit_length_in_cgs * const_unit_velocity_in_cgs
-l = l.to(l_units_cgs*unyt.g*unyt.cm**2/unyt.s).value
-L_r =np.array([ np.sum(l[rp_kpc<=r_values[i]],axis=0) for i in range(len(r_values))])
-norm_L_r = np.sqrt(L_r[:,0]**2+L_r[:,1]**2+L_r[:,2]**2)
-axs[3].scatter(r_values,norm_L_r,color='black', marker='+',s=100)
+l_units_cgs = (
+ const_unit_mass_in_cgs * const_unit_length_in_cgs * const_unit_velocity_in_cgs
+)
+l = l.to(l_units_cgs * unyt.g * unyt.cm ** 2 / unyt.s).value
+L_r = np.array([np.sum(l[rp_kpc <= r_values[i]], axis=0) for i in range(len(r_values))])
+norm_L_r = np.sqrt(L_r[:, 0] ** 2 + L_r[:, 1] ** 2 + L_r[:, 2] ** 2)
+# axs[3].scatter(r_values,norm_L_r,color='black', marker='+',s=100)
# NFW-like gas density profile
-def rho_r(r_value,f_b,M_200_cgs,r_200_cgs,c_200):
- rho_0 = M_200_cgs/(np.log(1+c_200)-c_200/(1+c_200))/(4*np.pi*r_200_cgs**3/c_200**3)
- result_cgs = rho_0*f_b/(c_200*r_value*(1+c_200*r_value)**2)
+def rho_r(r_value, f_b, M_200_cgs, r_200_cgs, c_200):
+ rho_0 = (
+ M_200_cgs
+ / (np.log(1 + c_200) - c_200 / (1 + c_200))
+ / (4 * np.pi * r_200_cgs ** 3 / c_200 ** 3)
+ )
+ result_cgs = rho_0 * f_b / (c_200 * r_value * (1 + c_200 * r_value) ** 2)
# Apply density cut
- rho_max_cgs = nH_max_cgs*m_H_cgs
+ rho_max_cgs = nH_max_cgs * m_H_cgs
result_cgs = np.array(result_cgs)
- result_cgs[result_cgs>rho_max_cgs]=rho_max_cgs
- return result_cgs
+ result_cgs[result_cgs > rho_max_cgs] = rho_max_cgs
+ return result_cgs
+
+
+# NFW-like gas mass inside a sphere with radius R
+def Mgas_r(r_value, f_b, M_200_cgs, r_200_cgs, c_200):
+ M_0 = M_200_cgs / (np.log(1 + c_200) - c_200 / (1 + c_200))
+ return (
+ M_0
+ * f_b
+ * (np.log(1 + c_200 * r_value) - c_200 * r_value / (1 + c_200 * r_value))
+ )
-# NFW-like gas mass inside a sphere with radius R
-def Mgas_r(r_value,f_b,M_200_cgs,r_200_cgs,c_200):
- M_0 = M_200_cgs/(np.log(1+c_200)-c_200/(1+c_200))
- return M_0 * f_b * (np.log(1+c_200*r_value)-c_200*r_value/(1+c_200*r_value))
# NFW Gravitational acceleration
-def a_NFW(r_value, M_200_cgs,r_200_cgs, c_200):
- a_pref = CONST_G_CGS*M_200_cgs/(np.log(1+c_200)-c_200/(1+c_200))/r_200_cgs**2
- return a_pref*((r_value/(r_value+1/c_200))-np.log(1+c_200*r_value))/r_value**2
+def a_NFW(r_value, M_200_cgs, r_200_cgs, c_200):
+ a_pref = (
+ CONST_G_CGS
+ * M_200_cgs
+ / (np.log(1 + c_200) - c_200 / (1 + c_200))
+ / r_200_cgs ** 2
+ )
+ return (
+ a_pref
+ * ((r_value / (r_value + 1 / c_200)) - np.log(1 + c_200 * r_value))
+ / r_value ** 2
+ )
+
# Integrate rho_gas*a_NFW
-def integrate(r_min, r_max, f_b, M_200_cgs, r_200_cgs, c_200, Nsteps = 10000):
+def integrate(r_min, r_max, f_b, M_200_cgs, r_200_cgs, c_200, Nsteps=10000):
# Perform the integration
r_range = np.linspace(r_min, r_max, Nsteps)
- dr = np.abs((r_max-r_min)/Nsteps)
- integrands = rho_r(r_range, f_b, M_200_cgs, r_200_cgs, c_200) * a_NFW(r_range, M_200_cgs,r_200_cgs, c_200)
- result_cgs = np.sum(integrands*dr)*r_200_cgs
+ dr = np.abs((r_max - r_min) / Nsteps)
+ integrands = rho_r(r_range, f_b, M_200_cgs, r_200_cgs, c_200) * a_NFW(
+ r_range, M_200_cgs, r_200_cgs, c_200
+ )
+ result_cgs = np.sum(integrands * dr) * r_200_cgs
return result_cgs
+
# NFW-like gas hydrostatic equilibrium internal energy profile
def u_vs_r(P0_cgs, r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200):
- result_cgs = (P0_cgs-integrate(r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200))/(gamma-1)/rho_r(r_value,f_b,M_200_cgs,r_200_cgs,c_200)
+ result_cgs = (
+ (P0_cgs - integrate(r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200))
+ / (gamma - 1)
+ / rho_r(r_value, f_b, M_200_cgs, r_200_cgs, c_200)
+ )
return result_cgs
+
# NFW-like gas hydrostatic equilibrium pressure profile
def P_vs_r(P0_cgs, r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200):
- result_cgs = P0_cgs-integrate(r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200)
+ result_cgs = P0_cgs - integrate(r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200)
return result_cgs
+
# NFW-like gas hydrostatic equilibrium temperature profile
-def T_vs_r(P0_cgs, r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200):
- result_cgs = u_vs_r(P0_cgs, r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200) * (gamma-1) * m_H_cgs/kb_cgs
+def T_vs_r(P0_cgs, r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200):
+ result_cgs = (
+ u_vs_r(P0_cgs, r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200)
+ * (gamma - 1)
+ * m_H_cgs
+ / kb_cgs
+ )
return result_cgs
-r_200_kpc = r_200_cgs / (PARSEC_IN_CGS*1e3)
-n_H_analytic = rho_r(np.abs(x)/r_200_kpc,f_b,M_200_cgs,r_200_cgs,c_200)/m_H_cgs
-M_gas_analytic = Mgas_r(np.abs(x)/r_200_kpc,f_b,M_200_cgs,r_200_cgs,c_200)/MSOL_IN_CGS
-# Angular momentum
+r_200_kpc = r_200_cgs / (PARSEC_IN_CGS * 1e3)
+n_H_analytic = rho_r(np.abs(x) / r_200_kpc, f_b, M_200_cgs, r_200_cgs, c_200) / m_H_cgs
+M_gas_analytic = (
+ Mgas_r(np.abs(x) / r_200_kpc, f_b, M_200_cgs, r_200_cgs, c_200) / MSOL_IN_CGS
+)
+
+# Angular momentum
Total_E = v_200 ** 2 / 2.0
-J_cgs = spin_lambda * const_G / np.sqrt(Total_E) * const_unit_length_in_cgs * const_unit_velocity_in_cgs
+J_cgs = (
+ spin_lambda
+ * const_G
+ / np.sqrt(Total_E)
+ * const_unit_length_in_cgs
+ * const_unit_velocity_in_cgs
+)
-L_200_cgs = J_cgs * M_200_cgs * f_b #np.linalg.norm(np.sum(l[rp_kpc<=r_200_kpc],axis=0))
-Mass_ratio = Mgas_r(np.abs(x)/r_200_kpc,f_b,M_200_cgs,r_200_cgs,c_200)/(f_b*M_200_cgs)
-L_analytic = Mass_ratio*L_200_cgs/l_units_cgs
+L_200_cgs = (
+ J_cgs * M_200_cgs * f_b
+) # np.linalg.norm(np.sum(l[rp_kpc<=r_200_kpc],axis=0))
+Mass_ratio = Mgas_r(np.abs(x) / r_200_kpc, f_b, M_200_cgs, r_200_cgs, c_200) / (
+ f_b * M_200_cgs
+)
+L_analytic = Mass_ratio * L_200_cgs / l_units_cgs
# Gas parameters
gamma = 5.0 / 3.0
-T0_cgs = T_200_cgs #1e5 # gas temperature on the edge of the box (if we want to set this manually)
-
-rho0_cgs = rho_r(Lbox_kpc/r_200_kpc*np.sqrt(3)/2,f_b,M_200_cgs,r_200_cgs,c_200) #gas density on the edge
-P0_cgs = rho0_cgs*kb_cgs*T0_cgs/m_H_cgs # gas pressure on the edge of the box
-
-Pressure_UNIT = MSOL_IN_CGS / PARSEC_IN_CGS / GYR_IN_CGS**2
-P_analytic = [P_vs_r(P0_cgs, np.abs(x[i]/r_200_kpc), Lbox_kpc/r_200_kpc*np.sqrt(3)/2, f_b, M_200_cgs, r_200_cgs, c_200)/Pressure_UNIT for i in range(map_pixel_length)] # gas particle internal energies
-#T_analytic = [T_vs_r(P0_cgs, np.abs(x[i]/r_200_kpc), Lbox_kpc/r_200_kpc*np.sqrt(3)/2, f_b, M_200_cgs, r_200_cgs, c_200) for i in range(map_pixel_length)] # gas particle internal energies
-
-
-#locs = [map_pixel_length / 4, map_pixel_length / 2, 3 * map_pixel_length / 4, map_pixel_length]
-#labels = np.array([Lbox_kpc / 4, Lbox_kpc / 2, 3*Lbox_kpc / 4, Lbox_kpc])
-locs = np.linspace(0,400,9)
-labels = np.round(locs,0)
+T0_cgs = (
+ T_200_cgs
+) # 1e5 # gas temperature on the edge of the box (if we want to set this manually)
+
+rho0_cgs = rho_r(
+ Lbox_kpc / r_200_kpc * np.sqrt(3) / 2, f_b, M_200_cgs, r_200_cgs, c_200
+) # gas density on the edge
+P0_cgs = rho0_cgs * kb_cgs * T0_cgs / m_H_cgs # gas pressure on the edge of the box
+
+Pressure_UNIT = MSOL_IN_CGS / PARSEC_IN_CGS / GYR_IN_CGS ** 2
+P_analytic = [
+ P_vs_r(
+ P0_cgs,
+ np.abs(x[i] / r_200_kpc),
+ Lbox_kpc / r_200_kpc * np.sqrt(3) / 2,
+ f_b,
+ M_200_cgs,
+ r_200_cgs,
+ c_200,
+ )
+ / Pressure_UNIT
+ for i in range(map_pixel_length)
+] # gas particle internal energies
+# T_analytic = [T_vs_r(P0_cgs, np.abs(x[i]/r_200_kpc), Lbox_kpc/r_200_kpc*np.sqrt(3)/2, f_b, M_200_cgs, r_200_cgs, c_200) for i in range(map_pixel_length)] # gas particle internal energies
+
+
+# locs = [map_pixel_length / 4, map_pixel_length / 2, 3 * map_pixel_length / 4, map_pixel_length]
+# labels = np.array([Lbox_kpc / 4, Lbox_kpc / 2, 3*Lbox_kpc / 4, Lbox_kpc])
+locs = np.linspace(0, 400, 9)
+labels = np.round(locs, 0)
axs[0].set_xlabel(r"$r$ [kpc]")
axs[0].set_xticks(locs, labels)
-axs[0].set_xlim(0, map_pixel_length/2)
-axs[0].plot(x, n_H_analytic, "k-",color='red')
+axs[0].set_xlim(0, map_pixel_length / 2)
+axs[0].plot(x, n_H_analytic, "k-", color="red")
axs[0].set_ylabel(r"$n_H(r)$ $[cm^{-3}]$")
-axs[0].axvline(x=r_200_kpc,color='gray',ls='dashed',label = '$R_200$')
-axs[0].set_xscale('log')
-axs[0].set_xlim(0.5,300)
+axs[0].axvline(x=r_200_kpc, color="gray", ls="dashed", label="$R_200$")
+axs[0].set_xscale("log")
+axs[0].set_xlim(0.5, 300)
axs[1].set_xlabel(r"$r$ [kpc]")
axs[1].set_xticks(locs, labels)
-axs[1].set_xlim(0, map_pixel_length/2)
-axs[1].plot(x, P_analytic, "k-",color='red')
+axs[1].set_xlim(0, map_pixel_length / 2)
+axs[1].plot(x, P_analytic, "k-", color="red")
axs[1].set_ylabel(r"$P(r)$ $[M_{sol} \cdot pc^{-1} \cdot Gyr^{-2}]$")
-axs[1].set_yscale('log')
+axs[1].set_yscale("log")
axs[1].set_yticks(np.logspace(3, 8, 6))
axs[1].set_ylim(1e3, 1e8)
-axs[1].axvline(x=r_200_kpc,color='gray',ls='dashed',label = '$R_200$')
-axs[1].set_xscale('log')
-axs[1].set_xlim(0.5,300)
+axs[1].axvline(x=r_200_kpc, color="gray", ls="dashed", label="$R_200$")
+axs[1].set_xscale("log")
+axs[1].set_xlim(0.5, 300)
axs[2].set_xlabel(r"$r$ [kpc]")
axs[2].set_xticks(locs, labels)
-axs[2].set_xlim(0, map_pixel_length/2)
-axs[2].plot(x, M_gas_analytic, "k-",color='red')
+axs[2].set_xlim(0, map_pixel_length / 2)
+axs[2].plot(x, M_gas_analytic, "k-", color="red")
axs[2].set_ylabel(r"$M_{gas}(<r)$ $[M_{sol}]$")
-axs[2].set_yscale('log')
+axs[2].set_yscale("log")
axs[2].set_yticks(np.logspace(7, 12, 6))
axs[2].set_ylim(1e7, 1e12)
-axs[2].axvline(x=r_200_kpc,color='gray',ls='dashed',label = '$R_200$')
-axs[2].axhline(y=f_b*M_200_cgs/(MSOL_IN_CGS),color='gray',ls='dashed',label = '$M_200$')
-axs[2].set_xscale('log')
-axs[2].set_xlim(0.5,300)
-
-axs[3].set_xlabel(r"$r$ [kpc]")
-axs[3].set_xticks(locs, labels)
-axs[3].set_xlim(0, map_pixel_length/2)
-axs[3].plot(x, L_analytic, "k-",color='red')
-axs[3].set_ylabel(r"$|\vec{L}_{gas}(<r)|$ $[M_{200}\cdot R_{200} \cdot v_{200} ]$ ")
-axs[3].set_yscale('log')
-axs[3].set_yticks(np.logspace(-4, 1, 6))
-axs[3].set_ylim(1e-4, 1e1)
-axs[3].axvline(x=r_200_kpc,color='gray',ls='dashed',label = '$R_200$')
-axs[3].set_xscale('log')
-axs[3].set_xlim(0.5,300)
-
-if with_reference:
- axs[0].plot(rho_vs_x_data[0],rho_vs_x_data[1],label='MFM $256^2$',color='red')
- #axs[1].plot(Bx_vs_x_data[0],Bx_vs_x_data[1],label='MFM $256^2$',color='red')
- #axs[2].plot(By_vs_x_data[0],By_vs_x_data[1],label='MFM $256^2$',color='red')
- #axs[3].plot(Berr_vs_x_data[0],Berr_vs_x_data[1],label='MFM $256^2$',color='red')
-
-#axs[2].legend()
-# Add panel with infromation about the run
-Ninfo = 4
+axs[2].axvline(x=r_200_kpc, color="gray", ls="dashed", label="$R_200$")
+axs[2].axhline(
+ y=f_b * M_200_cgs / (MSOL_IN_CGS), color="gray", ls="dashed", label="$M_200$"
+)
+axs[2].set_xscale("log")
+axs[2].set_xlim(0.5, 300)
+
+
+Ninfo = 3
text_common_args = dict(
fontsize=10, ha="center", va="center", transform=axs[Ninfo].transAxes
)
@@ -364,7 +407,6 @@ axs[Ninfo].text(
axs[Ninfo].axis("off")
-
for ax in axs:
ax.minorticks_on()
ax.grid()
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/README b/examples/MHDTests/CoolingHaloWithSpin_MHD/README
index 2c61dadafece03272185787a7ab6d148f13a6ddb..adc8ad69e6a50d981c11eb868785b30a1be9e552 100644
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/README
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/README
@@ -32,7 +32,7 @@ for runs with lambda cooling
Setup
-----
-There are two setups present: 1) dwarfsize, one similar to Peng Wang, Tom Abel, 0712.0872; 2) MWsize, similar to R. Pakmor, V. Springel, 1212.1452
+We follow R. Pakmor, V. Springel, 1212.1452
Snapshot slices
---------------
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/cooling_halo_MWsize.yml b/examples/MHDTests/CoolingHaloWithSpin_MHD/cooling_halo.yml
similarity index 100%
rename from examples/MHDTests/CoolingHaloWithSpin_MHD/cooling_halo_MWsize.yml
rename to examples/MHDTests/CoolingHaloWithSpin_MHD/cooling_halo.yml
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/cooling_halo_dwarfsize.yml b/examples/MHDTests/CoolingHaloWithSpin_MHD/cooling_halo_dwarfsize.yml
deleted file mode 100644
index aa7ce339bb875729b2a48b860b626daf81a227a3..0000000000000000000000000000000000000000
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/cooling_halo_dwarfsize.yml
+++ /dev/null
@@ -1,108 +0,0 @@
-# Define the system of units to use internally.
-InternalUnitSystem:
- UnitMass_in_cgs: 1.98848e39 # 10^6 solar masses
- UnitLength_in_cgs: 3.0856776e21 # Kiloparsecs
- UnitVelocity_in_cgs: 1e5 # Kilometres per second
- UnitCurrent_in_cgs: 2.088431e13 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
-
-# Parameters governing the time integration
-TimeIntegration:
- time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 10. # The end time of the simulation (in internal units).
- dt_min: 1e-12 #1e-8 # The minimal time-step size of the simulation (in internal units).
- dt_max: 1e-1 # The maximal time-step size of the simulation (in internal units).
-
-# Parameters governing the conserved quantities statistics
-Statistics:
- delta_time: 1e-2 # Time between statistics output
-
-Scheduler:
- tasks_per_cell: 200
-
-# Parameters governing the snapshots
-Snapshots:
- basename: CoolingHalo # Common part of the name of output files
- time_first: 0. # Time of the first output (in internal units)
- delta_time: 0.1 #0.1 # Time difference between consecutive outputs (in internal units)
-
-# Parameters for the hydrodynamics scheme
-SPH:
- resolution_eta: 1.3663 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
- minimal_temperature: 1e2 # Kelvin
-
-# Parameters for MHD schemes
-MHD:
- monopole_subtraction: 1.0
- artificial_diffusion: 0.0
- hyperbolic_dedner: 0.1
- hyperbolic_dedner_divv: 0.0
- parabolic_dedner: 1.0
- resistive_eta: 0.0
-
-# Parameters related to the initial conditions
-InitialConditions:
- file_name: CoolingHalo.hdf5 # The file to read
- periodic: 1
-
-# External potential parameters
-NFWPotential:
- useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
- position: [0.,0.,0.] # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
- M_200: 10000 # M200 of the galaxy disk
- h: 0.704 # reduced Hubble constant (value does not specify the used units!)
- concentration: 10.0 # concentration of the Halo
- timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
- epsilon: 0.01 # Softening size (internal units)
- #diskfraction: 0.0 # Disk mass fraction
- #bulgefraction: 0.0 # Bulge mass fraction
-
-# Parameters for the self-gravity scheme
-Gravity:
- eta: 0.025
- MAC: adaptive
- theta_cr: 0.7
- epsilon_fmm: 0.001
- max_physical_baryon_softening: 0.7 # Physical softening length (in internal units).
- mesh_side_length: 64
-
-# Simple model cooling parameters
-LambdaCooling:
- lambda_nH2_cgs: 1e-23 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3])
- cooling_tstep_mult: 0.1 # Dimensionless pre-factor for the time-step condition
-
-# Parameters for the EAGLE "equation of state"
-EAGLEEntropyFloor:
- Jeans_density_threshold_H_p_cm3: 1e-4 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
- Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
- Jeans_temperature_norm_K: 800 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
- Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
- Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
- Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
- Cool_temperature_norm_K: 10. # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin. (NOTE: This is below the min T and hence this floor does nothing)
- Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
-
-
-# Use solar abundances
-EAGLEChemistry:
- init_abundance_metal: 0.0129
- init_abundance_Hydrogen: 0.7065
- init_abundance_Helium: 0.2806
- init_abundance_Carbon: 0.00207
- init_abundance_Nitrogen: 0.000836
- init_abundance_Oxygen: 0.00549
- init_abundance_Neon: 0.00141
- init_abundance_Magnesium: 0.000591
- init_abundance_Silicon: 0.000683
- init_abundance_Iron: 0.0011
-
-# Standard EAGLE cooling options
-EAGLECooling:
- dir_name: ./coolingtables/ # Location of the Wiersma+09 cooling tables
- H_reion_z: 7.5 # Redshift of Hydrogen re-ionization
- H_reion_eV_p_H: 2.0 # Energy inject by Hydrogen re-ionization in electron-volt per Hydrogen atom
- He_reion_z_centre: 3.5 # Redshift of the centre of the Helium re-ionization Gaussian
- He_reion_z_sigma: 0.5 # Spread in redshift of the Helium re-ionization Gaussian
- He_reion_eV_p_H: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
-
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC.py b/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC.py
index c791b6c9c080f1d3d6be63579c147b727b320ffd..ecd957484ca3440b16b2d42c9add49b0025e5d56 100644
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC.py
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC.py
@@ -33,45 +33,55 @@ PARSEC_IN_CGS = 3.0856776e18
KM_PER_SEC_IN_CGS = 1.0e5
CONST_G_CGS = 6.672e-8
CONST_MU0_CGS = 4 * np.pi * 1e-2
-MSOL_IN_CGS = 1.9891e33 # Solar mass
-kb_cgs = 1.38e-16 # boltzmann constant
-m_H_cgs = 1.68e-24 # atomic hydrogen mass
+MSOL_IN_CGS = 1.9891e33 # Solar mass
+kb_cgs = 1.38e-16 # boltzmann constant
+m_H_cgs = 1.68e-24 # atomic hydrogen mass
# First set unit velocity and then the circular velocity parameter for the isothermal potential
const_unit_velocity_in_cgs = 1.0e5 # kms^-1
# Cosmological parameters
OMEGA = 0.3 # Cosmological matter fraction at z = 0
-h = 0.681 #0.67777 # hubble parameter
+h = 0.681 # 0.67777 # hubble parameter
# Find H_0, the inverse Hubble time, in cgs
H_0_cgs = 100.0 * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
# DM halo parameters
spin_lambda = 0.05 # spin parameter
+spin_lambda_choice = "Bullock" # which definition of spin parameter to use
f_b = 0.17 # baryon fraction
-c_200 = 7.2 # concentration parameter
+c_200 = 7.2 # concentration parameter
# Set the magnitude of the uniform seed magnetic field
-B0_Gaussian_Units = 1e-9 #1e-6 # 1 micro Gauss
+B0_Gaussian_Units = 1e-9 # 1e-6 # 1 micro Gauss
B0_cgs = np.sqrt(CONST_MU0_CGS / (4.0 * np.pi)) * B0_Gaussian_Units
# SPH
-eta = 1.3663 # kernel smoothing
+eta = 1.3663 # kernel smoothing
+
+# Additional options
+spin_lambda_choice = (
+ "Bullock"
+) # which definition of spin parameter to use (Peebles or Bullock)
+save_to_vtk = False
+plot_v_distribution = False
# From this we can find the virial radius, the radius within which the average density of the halo is
# 200. * the mean matter density
# Set M200 and get R200 and V200
-M_200_cgs = 1e12 * MSOL_IN_CGS
-rhoc_cgs = 3*H_0_cgs**2/(8*np.pi*CONST_G_CGS)
-r_200_cgs = (3*M_200_cgs/(4*np.pi*rhoc_cgs*200))**(1/3)
-v_200_cgs = np.sqrt(CONST_G_CGS*M_200_cgs/r_200_cgs)
-v_200 = v_200_cgs / const_unit_velocity_in_cgs
-T_200_cgs = m_H_cgs*v_200_cgs**2/(2*kb_cgs)
+M_200_cgs = 1e12 * MSOL_IN_CGS
+rhoc_cgs = 3 * H_0_cgs ** 2 / (8 * np.pi * CONST_G_CGS)
+r_200_cgs = (3 * M_200_cgs / (4 * np.pi * rhoc_cgs * 200)) ** (1 / 3)
+v_200_cgs = np.sqrt(CONST_G_CGS * M_200_cgs / r_200_cgs)
+v_200 = v_200_cgs / const_unit_velocity_in_cgs
+T_200_cgs = m_H_cgs * v_200_cgs ** 2 / (2 * kb_cgs)
# Gas parameters
gamma = 5.0 / 3.0
-T0_cgs = T_200_cgs #1e5 # gas temperature on the edge of the box (if we want to set this manually)
-nH_max_cgs = 1e0 # maximum hydrogen number density
+T0_cgs = (
+ T_200_cgs
+) # 1e5 # gas temperature on the edge of the box (if we want to set this manually)
+nH_max_cgs = 1e0 # maximum hydrogen number density
print("T_200 = %E" % T_200_cgs)
# Now set the unit length and mass
@@ -105,7 +115,7 @@ periodic = 1 # 1 For periodic box
boxSize = 4.0
G = const_G
N = int(sys.argv[1]) # Number of particles
-IAsource = 'IAfile'#'grid' #'IAfile'
+IAsource = "IAfile"
# Create the file
filename = "CoolingHalo.hdf5"
@@ -123,7 +133,7 @@ grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
# Loading initial arrangement file
-IAfile = 'glassCube_64.hdf5'
+IAfile = "glassCube_64.hdf5"
# smoothing kernel optimized for numpy Price 1012.1885 (6)
def open_IAfile(path_to_file):
IAfile = h5py.File(path_to_file, "r")
@@ -131,144 +141,165 @@ def open_IAfile(path_to_file):
h = IAfile["/PartType0/SmoothingLength"][:]
return pos, h
-if IAsource == 'IAfile':
+
+if IAsource == "IAfile":
# Loading initial arrangement file
- IAfile = 'glassCube_128.hdf5'
- coords,_ = open_IAfile(IAfile)
- lcut = (N/len(coords))**(1/3)
+ IAfile = "glassCube_128.hdf5"
+ coords, _ = open_IAfile(IAfile)
+ lcut = (N / len(coords)) ** (1 / 3)
coords -= 0.5
- mask = ((np.abs(coords[:,0])<=lcut/2) & (np.abs(coords[:,1])<=lcut/2) & (np.abs(coords[:,2])<=lcut/2))
+ mask = (
+ (np.abs(coords[:, 0]) <= lcut / 2)
+ & (np.abs(coords[:, 1]) <= lcut / 2)
+ & (np.abs(coords[:, 2]) <= lcut / 2)
+ )
coords = coords[mask]
coords /= lcut
- coords+=0.5
-
-elif IAsource == 'grid':
- Nside = int(N**(1/3))
+ coords += 0.5
+elif IAsource == "grid":
+ Nside = int(N ** (1 / 3))
grid_1d = np.linspace(0.0, 1.0, Nside) # Cube side is 1, centered at (0.5,0.5,0.5)
# Create a 3D grid of coordinates
x, y, z = np.meshgrid(grid_1d, grid_1d, grid_1d, indexing="ij")
# Reshape into a list of points
coords = np.vstack([x.ravel(), y.ravel(), z.ravel()]).T
- coords += 0.5/Nside
+ coords += 0.5 / Nside
+# functions for coordinate transformation
import scipy.optimize as opt
+
+
def nfw_cdf(r, r_s):
"""Computes the cumulative mass fraction of the NFW profile."""
- return (np.log(1 + r/r_s) - (r/r_s) / (1 + r/r_s))
+ return np.log(1 + r / r_s) - (r / r_s) / (1 + r / r_s)
-def invert_nfw_cdf(F_uni, r_s, R_max, Rmin = 1e-6):
+
+def invert_nfw_cdf(F_uni, r_s, R_max, Rmin=1e-6):
"""Find r_NFW by solving F_uni = F_NFW using a numerical solver."""
F_NFW_max = nfw_cdf(R_max, r_s) # Normalization factor
+
def objective(r_nfw):
return nfw_cdf(r_nfw, r_s) / F_NFW_max - F_uni # Find r where CDFs match
- return opt.root_scalar(objective, bracket=[Rmin, R_max], method='bisect').root
+
+ return opt.root_scalar(objective, bracket=[Rmin, R_max], method="bisect").root
+
# make center at zero
coords -= 0.5
# cut a sphere
r_uni = np.linalg.norm(coords, axis=1)
-coords = coords[r_uni<=0.5]
-
+coords = coords[r_uni <= 0.5]
coords *= boxSize * np.sqrt(3)
-# Save particle arrangement to .vtk file to view
-import vtk
-vtk_points = vtk.vtkPoints()
-for x, y, z in coords:
- vtk_points.InsertNextPoint(x, y, z)
-
-poly_data = vtk.vtkPolyData()
-poly_data.SetPoints(vtk_points)
-
-writer = vtk.vtkPolyDataWriter()
-writer.SetFileName("input.vtk")
-writer.SetInputData(poly_data)
-writer.Write() # you can use paraview to watch how particles are arranged
-
# calculate max distance from the center (units of r_200)
-R_max = np.sqrt(3)*(boxSize/2)
-r_s = 1/c_200
+R_max = np.sqrt(3) * (boxSize / 2)
+r_s = 1 / c_200
# calculate distances to the center
r_uni = np.linalg.norm(coords, axis=1)
# Compute uniform CDF
-F_uni = r_uni**3 / R_max**3
+F_uni = r_uni ** 3 / R_max ** 3
# Compute corresponding NFW radii
r_nfw = np.array([invert_nfw_cdf(F, r_s, R_max) for F in F_uni])
# Rescale positions
scaling_factors = r_nfw / r_uni
coords *= scaling_factors[:, np.newaxis]
-# Save particle arrangement to .vtk file to view
-import vtk
-vtk_points = vtk.vtkPoints()
-for x, y, z in coords:
- vtk_points.InsertNextPoint(x, y, z)
-
-poly_data = vtk.vtkPolyData()
-poly_data.SetPoints(vtk_points)
-
-writer = vtk.vtkPolyDataWriter()
-writer.SetFileName("rescaled.vtk")
-writer.SetInputData(poly_data)
-writer.Write() # you can use paraview to watch how particles are arranged
-
-
-
# NFW-like gas density profile
-def rho_r(r_value,f_b,M_200_cgs,r_200_cgs,c_200):
- rho_0 = M_200_cgs/(np.log(1+c_200)-c_200/(1+c_200))/(4*np.pi*r_200_cgs**3/c_200**3)
- result_cgs = rho_0*f_b/(c_200*r_value*(1+c_200*r_value)**2)
+def rho_r(r_value, f_b, M_200_cgs, r_200_cgs, c_200):
+ rho_0 = (
+ M_200_cgs
+ / (np.log(1 + c_200) - c_200 / (1 + c_200))
+ / (4 * np.pi * r_200_cgs ** 3 / c_200 ** 3)
+ )
+ result_cgs = rho_0 * f_b / (c_200 * r_value * (1 + c_200 * r_value) ** 2)
# Apply density cut
- rho_max_cgs = nH_max_cgs*m_H_cgs
+ rho_max_cgs = nH_max_cgs * m_H_cgs
result_cgs = np.array(result_cgs)
- result_cgs[result_cgs>rho_max_cgs]=rho_max_cgs
- return result_cgs
+ result_cgs[result_cgs > rho_max_cgs] = rho_max_cgs
+ return result_cgs
+
+
+# NFW-like gas mass inside a sphere with radius R
+def Mgas_r(r_value, f_b, M_200_cgs, r_200_cgs, c_200):
+ M_0 = M_200_cgs / (np.log(1 + c_200) - c_200 / (1 + c_200))
+ return (
+ M_0
+ * f_b
+ * (np.log(1 + c_200 * r_value) - c_200 * r_value / (1 + c_200 * r_value))
+ )
-# NFW-like gas mass inside a sphere with radius R
-def Mgas_r(r_value,f_b,M_200_cgs,r_200_cgs,c_200):
- M_0 = M_200_cgs/(np.log(1+c_200)-c_200/(1+c_200))
- return M_0 * f_b * (np.log(1+c_200*r_value)-c_200*r_value/(1+c_200*r_value))
# NFW Gravitational acceleration
-def a_NFW(r_value, M_200_cgs,r_200_cgs, c_200):
- a_pref = CONST_G_CGS*M_200_cgs/(np.log(1+c_200)-c_200/(1+c_200))/r_200_cgs**2
- return a_pref*((r_value/(r_value+1/c_200))-np.log(1+c_200*r_value))/r_value**2
+def a_NFW(r_value, M_200_cgs, r_200_cgs, c_200):
+ a_pref = (
+ CONST_G_CGS
+ * M_200_cgs
+ / (np.log(1 + c_200) - c_200 / (1 + c_200))
+ / r_200_cgs ** 2
+ )
+ return (
+ a_pref
+ * ((r_value / (r_value + 1 / c_200)) - np.log(1 + c_200 * r_value))
+ / r_value ** 2
+ )
+
# NFW Gravitational potential
-def phi_NFW(r_value, M_200_cgs,r_200_cgs, c_200):
- phi_pref = CONST_G_CGS*M_200_cgs/(np.log(1+c_200)-c_200/(1+c_200))/r_200_cgs
- return -phi_pref*np.log(1+c_200*r_value)/r_value
+def phi_NFW(r_value, M_200_cgs, r_200_cgs, c_200):
+ phi_pref = (
+ CONST_G_CGS * M_200_cgs / (np.log(1 + c_200) - c_200 / (1 + c_200)) / r_200_cgs
+ )
+ return -phi_pref * np.log(1 + c_200 * r_value) / r_value
+
# Integrate rho_gas*a_NFW
-def integrate(r_min, r_max, f_b, M_200_cgs, r_200_cgs, c_200, Nsteps = 10000):
+def integrate(r_min, r_max, f_b, M_200_cgs, r_200_cgs, c_200, Nsteps=10000):
# Perform the integration
r_range = np.linspace(r_min, r_max, Nsteps)
- dr = np.abs((r_max-r_min)/Nsteps)
- integrands = rho_r(r_range, f_b, M_200_cgs, r_200_cgs, c_200) * a_NFW(r_range, M_200_cgs,r_200_cgs, c_200)
- result_cgs = np.sum(integrands*dr)*r_200_cgs
+ dr = np.abs((r_max - r_min) / Nsteps)
+ integrands = rho_r(r_range, f_b, M_200_cgs, r_200_cgs, c_200) * a_NFW(
+ r_range, M_200_cgs, r_200_cgs, c_200
+ )
+ result_cgs = np.sum(integrands * dr) * r_200_cgs
return result_cgs
+
# NFW-like gas hydrostatic equilibrium internal energy profile
def u_vs_r(P0_cgs, r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200):
- result_cgs = (P0_cgs-integrate(r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200))/(gamma-1)/rho_r(r_value,f_b,M_200_cgs,r_200_cgs,c_200)
+ result_cgs = (
+ (P0_cgs - integrate(r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200))
+ / (gamma - 1)
+ / rho_r(r_value, f_b, M_200_cgs, r_200_cgs, c_200)
+ )
return result_cgs
+
# NFW-like gas hydrostatic equilibrium temperature profile
-def T_vs_r(P0_cgs, r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200):
- result_cgs = u_vs_r(P0_cgs, r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200) * (gamma-1) * m_H_cgs/kb_cgs
+def T_vs_r(P0_cgs, r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200):
+ result_cgs = (
+ u_vs_r(P0_cgs, r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200)
+ * (gamma - 1)
+ * m_H_cgs
+ / kb_cgs
+ )
return result_cgs
+
N = len(coords)
gas_mass = (
- Mgas_r(R_max,f_b,M_200_cgs,r_200_cgs,c_200)/M_200_cgs
+ Mgas_r(R_max, f_b, M_200_cgs, r_200_cgs, c_200) / M_200_cgs
) # get total gas mass within a sphere of radius sqrt(3)/2*Lbox
gas_particle_mass = gas_mass / float(N)
-print('Gas particle mass is %E ' % (gas_particle_mass*M_200_cgs/MSOL_IN_CGS))
+print("Gas particle mass is %E " % (gas_particle_mass * M_200_cgs / MSOL_IN_CGS))
# Unnormalized mass shell distribution
-r = np.logspace(np.log10(1e-6*boxSize),np.log10(boxSize * np.sqrt(3.0) / 2.0),round(10*N**(1/3)))
+r = np.logspace(
+ np.log10(1e-6 * boxSize),
+ np.log10(boxSize * np.sqrt(3.0) / 2.0),
+ round(10 * N ** (1 / 3)),
+)
# shift to centre of box
coords += np.full(coords.shape, boxSize / 2.0)
@@ -329,20 +360,19 @@ coords[:, 1] = y_coords
coords[:, 2] = z_coords
-# Save particle arrangement to .vtk file to view
-import vtk
-vtk_points = vtk.vtkPoints()
-for x, y, z in (coords-boxSize/2):
- vtk_points.InsertNextPoint(x, y, z)
-
-poly_data = vtk.vtkPolyData()
-poly_data.SetPoints(vtk_points)
-
-writer = vtk.vtkPolyDataWriter()
-writer.SetFileName("output.vtk")
-writer.SetInputData(poly_data)
-writer.Write() # you can use paraview to watch how particles are arranged
+# Save particle arrangement to .vtk file to open with paraview
+if save_to_vtk:
+ import vtk
+ vtk_points = vtk.vtkPoints()
+ for x, y, z in coords - boxSize / 2:
+ vtk_points.InsertNextPoint(x, y, z)
+ poly_data = vtk.vtkPolyData()
+ poly_data.SetPoints(vtk_points)
+ writer = vtk.vtkPolyDataWriter()
+ writer.SetFileName("output.vtk")
+ writer.SetInputData(poly_data)
+ writer.Write() # one can use paraview to watch how particles are arranged
radius = np.sqrt(
(coords[:, 0] - boxSize / 2.0) ** 2
@@ -350,8 +380,7 @@ radius = np.sqrt(
+ (coords[:, 2] - boxSize / 2.0) ** 2
)
axis_distance = np.sqrt(
- (coords[:, 0] - boxSize / 2.0) ** 2
- + (coords[:, 1] - boxSize / 2.0) ** 2
+ (coords[:, 0] - boxSize / 2.0) ** 2 + (coords[:, 1] - boxSize / 2.0) ** 2
)
# now give particle's velocities and magnetic fields
@@ -361,56 +390,68 @@ B = np.zeros((N, 3))
# first work out total angular momentum of the halo within the virial radius
# we work in units where r_vir = 1 and M_vir = 1
Total_E = v_200 ** 2 / 2.0
-#J = spin_lambda * const_G / np.sqrt(Total_E)
+# J = spin_lambda * const_G / np.sqrt(Total_E)
j_sp = spin_lambda * np.sqrt(2) * v_200 * 1
-print("j_sp =",j_sp )
+print("j_sp =", j_sp)
# all particles within the virial radius have omega parallel to the z-axis, magnitude
-# is proportional to 1 over the radius
+# is proportional to j/r^2
-def j(r_value,j_max,s,f_b,M_200_cgs,r_200_cgs,c_200):
- return j_max * (Mgas_r(r_value,f_b,M_200_cgs,r_200_cgs,c_200)/(M_200_cgs*f_b))**s
+# define specific angular momentum distribution
+def j(r_value, j_max, s, f_b, M_200_cgs, r_200_cgs, c_200):
+ return (
+ j_max
+ * (Mgas_r(r_value, f_b, M_200_cgs, r_200_cgs, c_200) / (M_200_cgs * f_b)) ** s
+ )
-def v_circ_r(r_value,f_b,M_200_cgs,r_200_cgs,c_200):
- return spin_lambda * v_200 * (a_NFW(r_value, M_200_cgs,r_200_cgs, c_200)/a_NFW(1, M_200_cgs,r_200_cgs, c_200) * r_value / 1)**0.5
omega = np.zeros((N, 3))
omega_ort = np.zeros((N, 3))
radius_vect = np.zeros((N, 3))
-l = np.zeros((N,3))
+l = np.zeros((N, 3))
for i in range(N):
- radius_vect[i,:] = coords[i, :] - boxSize / 2.0
- omega[i, 2] = j(radius[i],1,1,f_b,M_200_cgs,r_200_cgs,c_200)/radius[i]**2
- v[i, :] = np.cross(omega[i, :], radius_vect[i,:])
- #v[i, :] /= np.linalg.norm(v[i, :])
- #v[i, :] *= v_circ_r(radius[i],f_b,M_200_cgs,r_200_cgs,c_200)
+ radius_vect[i, :] = coords[i, :] - boxSize / 2.0
+ omega[i, 2] = j(radius[i], 1, 1, f_b, M_200_cgs, r_200_cgs, c_200) / radius[i] ** 2
+ v[i, :] = np.cross(omega[i, :], radius_vect[i, :])
B[i, 0] = B0_cgs / const_unit_magnetic_field_in_cgs
- l[i, :] = gas_particle_mass * np.cross(radius_vect[i,:],v[i, :])
+ l[i, :] = gas_particle_mass * np.cross(radius_vect[i, :], v[i, :])
+# select particles inside R_200
mask_r_200 = radius <= 1
-
-Lp_tot = np.sum(l[mask_r_200],axis=0)
+Lp_tot = np.sum(l[mask_r_200], axis=0)
Lp_tot_abs = np.linalg.norm(Lp_tot)
-
-normV = np.linalg.norm(v,axis = 1)
-
-Ep_kin = gas_particle_mass * normV**2 / 2
-Ep_pot = gas_particle_mass * phi_NFW(radius, M_200_cgs,r_200_cgs, c_200)/(const_unit_velocity_in_cgs**2)
-Ep_tot = np.sum(Ep_kin[mask_r_200]+Ep_pot[mask_r_200])
-
-calc_spin_par_P = Lp_tot_abs * np.sqrt(np.abs(Ep_tot))/const_G / (f_b*1)**(5/2)
-jp_sp = Lp_tot_abs/(gas_particle_mass*np.sum(mask_r_200))
-calc_spin_par_B = jp_sp/(np.sqrt(2) * 1 * v_200)
-v *= spin_lambda/calc_spin_par_B
-l *=spin_lambda/calc_spin_par_B
-
-normV = np.linalg.norm(v,axis = 1)
-import matplotlib.pyplot as plt
-fig,ax = plt.subplots(figsize = (5,5))
-ax.scatter(radius*r_200_cgs/(1e3*PARSEC_IN_CGS), normV)
-#ax.scatter(radius*r_200_cgs/(1e3*PARSEC_IN_CGS), v_circ_r(radius,f_b,M_200_cgs,r_200_cgs,c_200)/spin_lambda)
-ax.set_xlim(0,40)
-ax.set_ylim(0,100)
-plt.savefig('v_distr.png')
+normV = np.linalg.norm(v, axis=1)
+
+if spin_lambda_choice == "Peebles":
+ # Normalize to Peebles spin parameter
+ Ep_kin = gas_particle_mass * normV ** 2 / 2
+ Ep_pot = (
+ gas_particle_mass
+ * phi_NFW(radius, M_200_cgs, r_200_cgs, c_200)
+ / (const_unit_velocity_in_cgs ** 2)
+ )
+ Ep_tot = np.sum(Ep_kin[mask_r_200] + Ep_pot[mask_r_200])
+ calc_spin_par_P = (
+ Lp_tot_abs * np.sqrt(np.abs(Ep_tot)) / const_G / (f_b * 1) ** (5 / 2)
+ )
+ v *= spin_lambda / calc_spin_par_P
+ l *= spin_lambda / calc_spin_par_P
+elif spin_lambda_choice == "Bullock":
+ # Normalize to Bullock
+ jp_sp = Lp_tot_abs / (gas_particle_mass * np.sum(mask_r_200))
+ calc_spin_par_B = jp_sp / (np.sqrt(2) * 1 * v_200)
+ v *= spin_lambda / calc_spin_par_B
+ l *= spin_lambda / calc_spin_par_B
+
+# Draw velocity distribution
+if plot_v_distribution:
+ import matplotlib.pyplot as plt
+
+ normV = np.linalg.norm(v, axis=1)
+ fig, ax = plt.subplots(figsize=(5, 5))
+ ax.scatter(radius * r_200_cgs / (1e3 * PARSEC_IN_CGS), normV)
+ ax.set_xlim(0, 40)
+ ax.set_ylim(0, 100)
+ plt.savefig("v_distribution.png")
# Header
grp = file.create_group("/Header")
@@ -455,11 +496,18 @@ ds = grp.create_dataset("SmoothingLength", (N,), "f")
ds[()] = h
h = np.zeros(1)
-# Internal energies
-rho0_cgs = rho_r(boxSize*np.sqrt(3)/2,f_b,M_200_cgs,r_200_cgs,c_200) #gas density on the edge
-P0_cgs = rho0_cgs*kb_cgs*T0_cgs/m_H_cgs # gas pressure on the edge of the box
-u = [u_vs_r(P0_cgs, radius[i], boxSize*np.sqrt(3)/2, f_b, M_200_cgs, r_200_cgs, c_200)/const_unit_velocity_in_cgs**2 for i in range(N)] # gas particle internal energies
-#u_vs_r(radius)
+# Internal energies in hydrostatic equilibrium
+rho0_cgs = rho_r(
+ boxSize * np.sqrt(3) / 2, f_b, M_200_cgs, r_200_cgs, c_200
+) # gas density on the edge
+P0_cgs = rho0_cgs * kb_cgs * T0_cgs / m_H_cgs # gas pressure on the edge of the box
+u = [
+ u_vs_r(
+ P0_cgs, radius[i], boxSize * np.sqrt(3) / 2, f_b, M_200_cgs, r_200_cgs, c_200
+ )
+ / const_unit_velocity_in_cgs ** 2
+ for i in range(N)
+] # gas particle internal energies
u = np.full((N,), u)
ds = grp.create_dataset("InternalEnergy", (N,), "f")
ds[()] = u
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC_MWsize.py b/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC_MWsize.py
deleted file mode 100644
index d9814622f53f40d623cf82241df016dc89f9e491..0000000000000000000000000000000000000000
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC_MWsize.py
+++ /dev/null
@@ -1,279 +0,0 @@
-###############################################################################
-# This file is part of SWIFT.
-# Copyright (c) 2016 Stefan Arridge (stefan.arridge@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/>.
-#
-##############################################################################
-
-# Halo parameters similar to R. Pakmor, V. Springel, 1212.1452
-
-import h5py
-import sys
-import numpy as np
-import math
-import random
-
-# Generates N particles in a spherically symmetric distribution with density profile ~r^(-2)
-# usage: python3 makeIC.py 1000: generate 1000 particles
-
-# Some constants
-
-OMEGA = 0.3 # Cosmological matter fraction at z = 0
-PARSEC_IN_CGS = 3.0856776e18
-KM_PER_SEC_IN_CGS = 1.0e5
-CONST_G_CGS = 6.672e-8
-CONST_MU0_CGS = 4 * np.pi * 1e-2
-MSOL_IN_CGS = 1.9891e33 # Solar mass
-h = 0.67777 # hubble parameter
-gamma = 5.0 / 3.0
-eta = 1.3663 # kernel smoothing
-spin_lambda = 0.05 # spin parameter
-f_b = 0.17 # baryon fraction
-
-# First set unit velocity and then the circular velocity parameter for the isothermal potential
-const_unit_velocity_in_cgs = 1.0e5 # kms^-1
-
-# Set the magnitude of the uniform seed magnetic field
-
-B0_Gaussian_Units = 1e-9 #1e-6 # 1 micro Gauss
-B0_cgs = np.sqrt(CONST_MU0_CGS / (4.0 * np.pi)) * B0_Gaussian_Units
-
-# Now we use this to get the virial mass and virial radius, which we will set to be the unit mass and radius
-
-# Find H_0, the inverse Hubble time, in cgs
-
-H_0_cgs = 100.0 * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
-
-# From this we can find the virial radius, the radius within which the average density of the halo is
-# 200. * the mean matter density
-
-
-# Set M200 and get R200 and V200
-M_200_cgs = 1e12 * MSOL_IN_CGS
-rhoc_cgs = 3*H_0_cgs**2/(8*np.pi*CONST_G_CGS)
-r_200_cgs = (3*M_200_cgs/(4*np.pi*rhoc_cgs))**(1/3)
-v_200_cgs = np.sqrt(CONST_G_CGS*M_200_cgs/r_200_cgs)
-v_200 = v_200_cgs / const_unit_velocity_in_cgs
-
-# Now set the unit length and mass
-
-const_unit_mass_in_cgs = M_200_cgs
-const_unit_length_in_cgs = r_200_cgs
-
-print("UnitMass_in_cgs: ", const_unit_mass_in_cgs)
-print("UnitLength_in_cgs: ", const_unit_length_in_cgs)
-print("UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs)
-
-# Now set the magnetic field unit
-
-const_unit_magnetic_field_in_cgs = 1e-7 # 1muG
-
-# Derived quantities
-const_unit_time_in_cgs = const_unit_length_in_cgs / const_unit_velocity_in_cgs
-const_unit_current_in_cgs = const_unit_mass_in_cgs / (
- const_unit_magnetic_field_in_cgs * const_unit_time_in_cgs ** 2
-)
-
-print("UnitTime_in_cgs: ", const_unit_time_in_cgs)
-print("UnitCurrent_in_cgs: ", const_unit_current_in_cgs)
-
-const_G = (
- CONST_G_CGS
- * const_unit_mass_in_cgs
- * const_unit_time_in_cgs
- * const_unit_time_in_cgs
- / (const_unit_length_in_cgs * const_unit_length_in_cgs * const_unit_length_in_cgs)
-)
-print("G=", const_G)
-
-# Parameters
-periodic = 1 # 1 For periodic box
-boxSize = 4.0
-G = const_G
-N = int(sys.argv[1]) # Number of particles
-
-# Create the file
-filename = "CoolingHalo.hdf5"
-file = h5py.File(filename, "w")
-
-# Units
-grp = file.create_group("/Units")
-grp.attrs["Unit length in cgs (U_L)"] = const_unit_length_in_cgs
-grp.attrs["Unit mass in cgs (U_M)"] = const_unit_mass_in_cgs
-grp.attrs["Unit time in cgs (U_t)"] = (
- const_unit_length_in_cgs / const_unit_velocity_in_cgs
-)
-grp.attrs["Unit current in cgs (U_I)"] = const_unit_current_in_cgs
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
-
-# set seed for random number
-np.random.seed(1234)
-
-gas_mass = (
- f_b * np.sqrt(3.0) / 2.0
-) # virial mass of halo is 1, virial radius is 1, enclosed mass scales with r
-gas_particle_mass = gas_mass / float(N)
-
-# Positions
-# r^(-2) distribution corresponds to uniform distribution in radius
-radius = (
- boxSize * np.sqrt(3.0) / 2.0 * np.random.rand(N)
-) # the diagonal extent of the cube
-ctheta = -1.0 + 2 * np.random.rand(N)
-stheta = np.sqrt(1.0 - ctheta ** 2)
-phi = 2 * math.pi * np.random.rand(N)
-coords = np.zeros((N, 3))
-coords[:, 0] = radius * stheta * np.cos(phi)
-coords[:, 1] = radius * stheta * np.sin(phi)
-coords[:, 2] = radius * ctheta
-
-# shift to centre of box
-coords += np.full((N, 3), boxSize / 2.0)
-print("x range = (%f,%f)" % (np.min(coords[:, 0]), np.max(coords[:, 0])))
-print("y range = (%f,%f)" % (np.min(coords[:, 1]), np.max(coords[:, 1])))
-print("z range = (%f,%f)" % (np.min(coords[:, 2]), np.max(coords[:, 2])))
-
-print(np.mean(coords[:, 0]))
-print(np.mean(coords[:, 1]))
-print(np.mean(coords[:, 2]))
-
-# now find the particles which are within the box
-
-x_coords = coords[:, 0]
-y_coords = coords[:, 1]
-z_coords = coords[:, 2]
-
-ind = np.where(x_coords < boxSize)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-ind = np.where(x_coords > 0.0)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-ind = np.where(y_coords < boxSize)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-ind = np.where(y_coords > 0.0)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-ind = np.where(z_coords < boxSize)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-ind = np.where(z_coords > 0.0)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-# count number of particles
-
-N = x_coords.size
-
-print("Number of particles in the box = ", N)
-
-# make the coords and radius arrays again
-coords = np.zeros((N, 3))
-coords[:, 0] = x_coords
-coords[:, 1] = y_coords
-coords[:, 2] = z_coords
-
-radius = np.sqrt(
- (coords[:, 0] - boxSize / 2.0) ** 2
- + (coords[:, 1] - boxSize / 2.0) ** 2
- + (coords[:, 2] - boxSize / 2.0) ** 2
-)
-
-# now give particle's velocities and magnetic fields
-v = np.zeros((N, 3))
-B = np.zeros((N, 3))
-
-# first work out total angular momentum of the halo within the virial radius
-# we work in units where r_vir = 1 and M_vir = 1
-Total_E = v_200 ** 2 / 2.0
-J = spin_lambda * const_G / np.sqrt(Total_E)
-print("J =", J)
-# all particles within the virial radius have omega parallel to the z-axis, magnitude
-# is proportional to 1 over the radius
-omega = np.zeros((N, 3))
-for i in range(N):
- omega[i, 2] = 3.0 * J / radius[i]
- v[i, :] = np.cross(omega[i, :], (coords[i, :] - boxSize / 2.0))
- #B[i, 2] = B0_cgs / const_unit_magnetic_field_in_cgs
- B[i, 0] = B0_cgs / const_unit_magnetic_field_in_cgs
-
-# Header
-grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = boxSize
-grp.attrs["NumPart_Total"] = [N, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_ThisFile"] = [N, 0, 0, 0, 0, 0]
-grp.attrs["Time"] = 0.0
-grp.attrs["NumFilesPerSnapshot"] = 1
-grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["Dimension"] = 3
-
-# Particle group
-grp = file.create_group("/PartType0")
-
-ds = grp.create_dataset("Coordinates", (N, 3), "d")
-ds[()] = coords
-coords = np.zeros(1)
-
-ds = grp.create_dataset("Velocities", (N, 3), "f")
-ds[()] = v
-v = np.zeros(1)
-
-# Magnetic field
-ds = grp.create_dataset("MagneticFluxDensities", (N, 3), "f")
-ds[()] = B
-B = np.zeros(1)
-
-# All particles of equal mass
-m = np.full((N,), gas_particle_mass)
-ds = grp.create_dataset("Masses", (N,), "f")
-ds[()] = m
-m = np.zeros(1)
-
-# Smoothing lengths
-l = (4.0 * np.pi * radius ** 2 / N) ** (
- 1.0 / 3.0
-) # local mean inter-particle separation
-h = np.full((N,), eta * l)
-ds = grp.create_dataset("SmoothingLength", (N,), "f")
-ds[()] = h
-h = np.zeros(1)
-
-# Internal energies
-u = v_200 ** 2 / (2.0 * (gamma - 1.0))
-u = np.full((N,), u)
-ds = grp.create_dataset("InternalEnergy", (N,), "f")
-ds[()] = u
-u = np.zeros(1)
-
-# Particle IDs
-ids = 1 + np.linspace(0, N, N, endpoint=False)
-ds = grp.create_dataset("ParticleIDs", (N,), "L")
-ds[()] = ids
-
-file.close()
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC_MWsize_NFWIC.py b/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC_MWsize_NFWIC.py
deleted file mode 100644
index 99ebb30a8824b161b1fcade44ec54d3e2525ef4c..0000000000000000000000000000000000000000
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC_MWsize_NFWIC.py
+++ /dev/null
@@ -1,411 +0,0 @@
-###############################################################################
-# This file is part of SWIFT.
-# Copyright (c) 2016 Stefan Arridge (stefan.arridge@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/>.
-#
-##############################################################################
-
-# Halo parameters similar to R. Pakmor, V. Springel, 1212.1452
-
-import h5py
-import sys
-import numpy as np
-import math
-import random
-
-# Generates N particles in a spherically symmetric distribution with density profile ~r^(-2)
-# usage: python3 makeIC.py 1000: generate 1000 particles
-
-# Some constants cgs
-PARSEC_IN_CGS = 3.0856776e18
-KM_PER_SEC_IN_CGS = 1.0e5
-CONST_G_CGS = 6.672e-8
-CONST_MU0_CGS = 4 * np.pi * 1e-2
-MSOL_IN_CGS = 1.9891e33 # Solar mass
-kb_cgs = 1.38e-16 # boltzmann constant
-m_H_cgs = 1.68e-24 # atomic hydrogen mass
-# First set unit velocity and then the circular velocity parameter for the isothermal potential
-const_unit_velocity_in_cgs = 1.0e5 # kms^-1
-
-# Cosmological parameters
-OMEGA = 0.3 # Cosmological matter fraction at z = 0
-h = 0.681 #0.67777 # hubble parameter
-# Find H_0, the inverse Hubble time, in cgs
-H_0_cgs = 100.0 * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
-
-# DM halo parameters
-spin_lambda = 0.05 # spin parameter
-f_b = 0.17 # baryon fraction
-c_200 = 7.2 # concentration parameter
-
-# Set the magnitude of the uniform seed magnetic field
-B0_Gaussian_Units = 1e-9 #1e-6 # 1 micro Gauss
-B0_cgs = np.sqrt(CONST_MU0_CGS / (4.0 * np.pi)) * B0_Gaussian_Units
-
-# SPH
-eta = 1.3663 # kernel smoothing
-
-# From this we can find the virial radius, the radius within which the average density of the halo is
-# 200. * the mean matter density
-
-# Set M200 and get R200 and V200
-M_200_cgs = 1e12 * MSOL_IN_CGS
-rhoc_cgs = 3*H_0_cgs**2/(8*np.pi*CONST_G_CGS)
-r_200_cgs = (3*M_200_cgs/(4*np.pi*rhoc_cgs*200))**(1/3)
-v_200_cgs = np.sqrt(CONST_G_CGS*M_200_cgs/r_200_cgs)
-v_200 = v_200_cgs / const_unit_velocity_in_cgs
-T_200_cgs = m_H_cgs*v_200_cgs**2/(2*kb_cgs)
-
-# Gas parameters
-gamma = 5.0 / 3.0
-T0_cgs = T_200_cgs #1e5 # gas temperature on the edge of the box (if we want to set this manually)
-nH_max_cgs = 1e0 # maximum hydrogen number density
-print("T_200 = %E" % T_200_cgs)
-
-# Now set the unit length and mass
-const_unit_mass_in_cgs = M_200_cgs
-const_unit_length_in_cgs = r_200_cgs
-print("UnitMass_in_cgs: ", const_unit_mass_in_cgs)
-print("UnitLength_in_cgs: ", const_unit_length_in_cgs)
-print("UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs)
-
-# Now set the magnetic field unit
-const_unit_magnetic_field_in_cgs = 1e-7 # 1muG
-
-# Derived quantities
-const_unit_time_in_cgs = const_unit_length_in_cgs / const_unit_velocity_in_cgs
-const_unit_current_in_cgs = const_unit_mass_in_cgs / (
- const_unit_magnetic_field_in_cgs * const_unit_time_in_cgs ** 2
-)
-print("UnitTime_in_cgs: ", const_unit_time_in_cgs)
-print("UnitCurrent_in_cgs: ", const_unit_current_in_cgs)
-const_G = (
- CONST_G_CGS
- * const_unit_mass_in_cgs
- * const_unit_time_in_cgs
- * const_unit_time_in_cgs
- / (const_unit_length_in_cgs * const_unit_length_in_cgs * const_unit_length_in_cgs)
-)
-print("G=", const_G)
-
-# Parameters
-periodic = 1 # 1 For periodic box
-boxSize = 4.0
-G = const_G
-N = int(sys.argv[1]) # Number of particles
-
-# Create the file
-filename = "CoolingHalo.hdf5"
-file = h5py.File(filename, "w")
-
-# Units
-grp = file.create_group("/Units")
-grp.attrs["Unit length in cgs (U_L)"] = const_unit_length_in_cgs
-grp.attrs["Unit mass in cgs (U_M)"] = const_unit_mass_in_cgs
-grp.attrs["Unit time in cgs (U_t)"] = (
- const_unit_length_in_cgs / const_unit_velocity_in_cgs
-)
-grp.attrs["Unit current in cgs (U_I)"] = const_unit_current_in_cgs
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
-
-# NFW-like gas density profile
-def rho_r(r_value,f_b,M_200_cgs,r_200_cgs,c_200):
- rho_0 = M_200_cgs/(np.log(1+c_200)-c_200/(1+c_200))/(4*np.pi*r_200_cgs**3/c_200**3)
- result_cgs = rho_0*f_b/(c_200*r_value*(1+c_200*r_value)**2)
- # Apply density cut
- rho_max_cgs = nH_max_cgs*m_H_cgs
- result_cgs = np.array(result_cgs)
- result_cgs[result_cgs>rho_max_cgs]=rho_max_cgs
- return result_cgs
-
-# NFW-like gas mass inside a sphere with radius R
-def Mgas_r(r_value,f_b,M_200_cgs,r_200_cgs,c_200):
- M_0 = M_200_cgs/(np.log(1+c_200)-c_200/(1+c_200))
- return M_0 * f_b * (np.log(1+c_200*r_value)-r_value/(1+c_200*r_value))
-
-# NFW Gravitational acceleration
-def a_NFW(r_value, M_200_cgs,r_200_cgs, c_200):
- a_pref = CONST_G_CGS*M_200_cgs/(np.log(1+c_200)-c_200/(1+c_200))/r_200_cgs**2
- return a_pref*((r_value/(r_value+1/c_200))-np.log(1+c_200*r_value))/r_value**2
-
-# Integrate rho_gas*a_NFW
-def integrate(r_min, r_max, f_b, M_200_cgs, r_200_cgs, c_200, Nsteps = 1000):
- # Perform the integration
- r_range = np.linspace(r_min, r_max, Nsteps)
- dr = np.abs((r_max-r_min)/Nsteps)
- integrands = rho_r(r_range, f_b, M_200_cgs, r_200_cgs, c_200) * a_NFW(r_range, M_200_cgs,r_200_cgs, c_200)
- result_cgs = np.sum(integrands*dr)*r_200_cgs
- return result_cgs
-
-# NFW-like gas hydrostatic equilibrium internal energy profile
-def u_vs_r(P0_cgs, r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200):
- result_cgs = (P0_cgs-integrate(r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200))/(gamma-1)/rho_r(r_value,f_b,M_200_cgs,r_200_cgs,c_200)
- return result_cgs
-
-
-# NFW-like gas hydrostatic equilibrium temperature profile
-def T_vs_r(P0_cgs, r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200):
- result_cgs = u_vs_r(P0_cgs, r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200) * (gamma-1) * m_H_cgs/kb_cgs
- return result_cgs
-
-# Masses
-gas_mass = (
- Mgas_r(boxSize*np.sqrt(3)/2,f_b,M_200_cgs,r_200_cgs,c_200)/M_200_cgs
-) # get total gas mass within a sphere of radius sqrt(3)/2*Lbox
-gas_particle_mass = gas_mass / float(N)
-print('Gas particle mass is %E ' % (gas_particle_mass*1e12))
-
-# Unnormalized mass shell distribution
-r = np.logspace(np.log10(1e-6*boxSize),np.log10(boxSize * np.sqrt(3.0) / 2.0),round(10*N**(1/3)))
-
-# fibonacci-based grid spherical distribution
-def fibonacci_sphere(Np):
- """Returns `Np` points distributed on a unit sphere using a Fibonacci lattice."""
- indicesp = np.arange(0, Np, dtype=float) + 0.5
- phip = (1 + np.sqrt(5)) / 2 # Golden ratio
- thetap = 2 * np.pi * indicesp / phip
- zp = 1 - (2 * indicesp / Np)
- mask_badzp = np.abs(zp)>1
- zp[mask_badzp] = np.sign(zp[mask_badzp])*1
- radiusp = np.sqrt(1 - zp**2)
- xp, yp = radiusp * np.cos(thetap), radiusp * np.sin(thetap)
-
- coordsp = np.vstack((xp, yp, zp)).T
- return coordsp
-
-# fill volume with spherical shells with a given dN vs r distribution
-def generate_coordinates(r_valuesp, N_vs_r_dr, noise_drp):
- """Generates 3D coordinates based on a given radial distribution and particle counts."""
- all_coords = []
-
- for i, (rp, Np, drp) in enumerate(zip(r_valuesp, N_vs_r_dr, noise_drp)):
- if Np > 2: # Avoid empty shells
- sphere_points = fibonacci_sphere(Np) # Generate points on unit sphere
- scaled_points = sphere_points * rp # Scale by radius
-
- # add noise to the shell
- noise_drp = np.random.randn(int(Np),3)
- noise_drp *= 1/np.linalg.norm(noise_drp, axis=1, keepdims=True) # Normalize
- noise_magnitudes = np.random.normal(loc=0, scale=drp/2, size=(int(Np), 1)) # Add gaussian noise amplitude distribution, with 2sigma = shell thickness
- noise_drp *= noise_magnitudes
- scaled_points += noise_drp
-
- # append shell
- all_coords.append(scaled_points)
-
- return np.vstack(all_coords) # Stack all coordinates into a single array
-
-# selects non-overlapping shells
-def minimal_covering_intervals(intervals, x_min, x_max):
-
- # Sort intervals based on their starting points (and ending points in case of ties)
-
- i = 0
- indeces = np.array([0])
- counter = 0
-
- while i<len(intervals)-1:
- distances_to_left = intervals[:,0]-intervals[i,1] # select ith interval and find distances to left edges
- indeces_non_intersect = np.argwhere(distances_to_left>0) # find intervals that don't intersect
- if len(indeces_non_intersect)>0: # if there are still non-intersecting intervals to the left
- ilast = np.min(indeces_non_intersect)-1 # find index of last intersecting interval
- else:
- break
-
- if i<ilast:
- indeces = np.append(indeces,ilast)
- i=ilast
- else:
- i=ilast+1
-
- # Loop safety
- counter+=1
- if counter>len(intervals):
- print('program loop')
- break
-
- return indeces
-
-# get particle distribution for the selected gas profile (NFW-like)
-def generate_particle_distribution(shell_sampling = round(10*N**(1/3)), noise_level = 1e-1):
-
- r = np.logspace(np.log10(1e-6*boxSize),np.log10(boxSize * np.sqrt(3.0) / 2.0),shell_sampling) # radial sampling
- d_bin = (gas_particle_mass*M_200_cgs/rho_r(r, f_b, M_200_cgs, r_200_cgs, c_200))**(1/3)/const_unit_length_in_cgs # estimate shell thickness
- dN_bin = np.round(4*np.pi*r**2/(d_bin)**2,0) # estimate amount of particles for each shell
- bins = np.vstack([r-d_bin/2,r+d_bin/2]).T # get shell radial intervals
- indeces_selected = minimal_covering_intervals(bins, d_bin[0]/2, r[-1]) # select non-overlapping shells
- coords = generate_coordinates(r[indeces_selected],dN_bin[indeces_selected],noise_level*d_bin[indeces_selected]) # generate particle distribution
-
- return coords
-
-coords = generate_particle_distribution(noise_level=0.25)
-
-# shift to centre of box
-coords += np.full(coords.shape, boxSize / 2.0)
-print("x range = (%f,%f)" % (np.min(coords[:, 0]), np.max(coords[:, 0])))
-print("y range = (%f,%f)" % (np.min(coords[:, 1]), np.max(coords[:, 1])))
-print("z range = (%f,%f)" % (np.min(coords[:, 2]), np.max(coords[:, 2])))
-
-print(np.mean(coords[:, 0]))
-print(np.mean(coords[:, 1]))
-print(np.mean(coords[:, 2]))
-
-# now find the particles which are within the box
-
-x_coords = coords[:, 0]
-y_coords = coords[:, 1]
-z_coords = coords[:, 2]
-
-ind = np.where(x_coords < boxSize)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-ind = np.where(x_coords > 0.0)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-ind = np.where(y_coords < boxSize)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-ind = np.where(y_coords > 0.0)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-ind = np.where(z_coords < boxSize)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-ind = np.where(z_coords > 0.0)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-# count number of particles
-
-N = x_coords.size
-
-print("Number of particles in the box = ", N)
-
-# make the coords and radius arrays again
-coords = np.zeros((N, 3))
-coords[:, 0] = x_coords
-coords[:, 1] = y_coords
-coords[:, 2] = z_coords
-
-
-# Save particle arrangement to .vtk file to view
-import vtk
-vtk_points = vtk.vtkPoints()
-for x, y, z in coords:
- vtk_points.InsertNextPoint(x, y, z)
-
-poly_data = vtk.vtkPolyData()
-poly_data.SetPoints(vtk_points)
-
-writer = vtk.vtkPolyDataWriter()
-writer.SetFileName("output.vtk")
-writer.SetInputData(poly_data)
-writer.Write() # you can use paraview to watch how particles are arranged
-
-
-radius = np.sqrt(
- (coords[:, 0] - boxSize / 2.0) ** 2
- + (coords[:, 1] - boxSize / 2.0) ** 2
- + (coords[:, 2] - boxSize / 2.0) ** 2
-)
-
-# now give particle's velocities and magnetic fields
-v = np.zeros((N, 3))
-B = np.zeros((N, 3))
-
-# first work out total angular momentum of the halo within the virial radius
-# we work in units where r_vir = 1 and M_vir = 1
-Total_E = v_200 ** 2 / 2.0
-J = spin_lambda * const_G / np.sqrt(Total_E)
-print("J =", J)
-# all particles within the virial radius have omega parallel to the z-axis, magnitude
-# is proportional to 1 over the radius
-omega = np.zeros((N, 3))
-for i in range(N):
- omega[i, 2] = 3.0 * J / radius[i]
- v[i, :] = np.cross(omega[i, :], (coords[i, :] - boxSize / 2.0))
- #B[i, 2] = B0_cgs / const_unit_magnetic_field_in_cgs
- B[i, 0] = B0_cgs / const_unit_magnetic_field_in_cgs
-
-# Header
-grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = boxSize
-grp.attrs["NumPart_Total"] = [N, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_ThisFile"] = [N, 0, 0, 0, 0, 0]
-grp.attrs["Time"] = 0.0
-grp.attrs["NumFilesPerSnapshot"] = 1
-grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["Dimension"] = 3
-
-# Particle group
-grp = file.create_group("/PartType0")
-
-ds = grp.create_dataset("Coordinates", (N, 3), "d")
-ds[()] = coords
-coords = np.zeros(1)
-
-ds = grp.create_dataset("Velocities", (N, 3), "f")
-ds[()] = v
-v = np.zeros(1)
-
-# Magnetic field
-ds = grp.create_dataset("MagneticFluxDensities", (N, 3), "f")
-ds[()] = B
-B = np.zeros(1)
-
-# All particles of equal mass
-m = np.full((N,), gas_particle_mass)
-ds = grp.create_dataset("Masses", (N,), "f")
-ds[()] = m
-m = np.zeros(1)
-
-# Smoothing lengths
-l = (4.0 * np.pi * radius ** 2 / N) ** (
- 1.0 / 3.0
-) # local mean inter-particle separation
-h = np.full((N,), eta * l)
-ds = grp.create_dataset("SmoothingLength", (N,), "f")
-ds[()] = h
-h = np.zeros(1)
-
-# Internal energies
-rho0_cgs = rho_r(boxSize*np.sqrt(3)/2,f_b,M_200_cgs,r_200_cgs,c_200) #gas density on the edge
-P0_cgs = rho0_cgs*kb_cgs*T0_cgs/m_H_cgs # gas pressure on the edge of the box
-u = [u_vs_r(P0_cgs, radius[i], boxSize*np.sqrt(3)/2, f_b, M_200_cgs, r_200_cgs, c_200)/const_unit_velocity_in_cgs**2 for i in range(N)] # gas particle internal energies
-#u_vs_r(radius)
-u = np.full((N,), u)
-ds = grp.create_dataset("InternalEnergy", (N,), "f")
-ds[()] = u
-u = np.zeros(1)
-
-# Particle IDs
-ids = 1 + np.linspace(0, N, N, endpoint=False)
-ds = grp.create_dataset("ParticleIDs", (N,), "L")
-ds[()] = ids
-
-file.close()
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC_dwarfsize.py b/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC_dwarfsize.py
deleted file mode 100644
index 299d3e09c525558f15cb1dcdb5115c84d0991d81..0000000000000000000000000000000000000000
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC_dwarfsize.py
+++ /dev/null
@@ -1,279 +0,0 @@
-###############################################################################
-# This file is part of SWIFT.
-# Copyright (c) 2016 Stefan Arridge (stefan.arridge@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/>.
-#
-##############################################################################
-
-# Halo parameters similar to Peng Wang, Tom Abel, 0712.0872
-
-import h5py
-import sys
-import numpy as np
-import math
-import random
-
-# Generates N particles in a spherically symmetric distribution with density profile ~r^(-2)
-# usage: python3 makeIC.py 1000: generate 1000 particles
-
-# Some constants
-
-OMEGA = 0.3 # Cosmological matter fraction at z = 0
-PARSEC_IN_CGS = 3.0856776e18
-KM_PER_SEC_IN_CGS = 1.0e5
-CONST_G_CGS = 6.672e-8
-CONST_MU0_CGS = 4 * np.pi * 1e-2
-MSOL_IN_CGS = 1.9891e33 # Solar mass
-h = 0.67777 # hubble parameter
-gamma = 5.0 / 3.0
-eta = 1.3663 # kernel smoothing
-spin_lambda = 0.05 # spin parameter
-f_b = 0.1 # baryon fraction
-
-# First set unit velocity and then the circular velocity parameter for the isothermal potential
-const_unit_velocity_in_cgs = 1.0e5 # kms^-1
-
-# Set the magnitude of the uniform seed magnetic field
-
-B0_Gaussian_Units = 1e-9 #1e-6 # 1 micro Gauss
-B0_cgs = np.sqrt(CONST_MU0_CGS / (4.0 * np.pi)) * B0_Gaussian_Units
-
-# Now we use this to get the virial mass and virial radius, which we will set to be the unit mass and radius
-
-# Find H_0, the inverse Hubble time, in cgs
-
-H_0_cgs = 100.0 * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
-
-# From this we can find the virial radius, the radius within which the average density of the halo is
-# 200. * the mean matter density
-
-
-# Set M200 and get R200 and V200
-M_200_cgs = 1e10 * MSOL_IN_CGS
-rhoc_cgs = 3*H_0_cgs**2/(8*np.pi*CONST_G_CGS)
-r_200_cgs = (3*M_200_cgs/(4*np.pi*rhoc_cgs))**(1/3)
-v_200_cgs = np.sqrt(CONST_G_CGS*M_200_cgs/r_200_cgs)
-v_200 = v_200_cgs / const_unit_velocity_in_cgs
-
-# Now set the unit length and mass
-
-const_unit_mass_in_cgs = M_200_cgs
-const_unit_length_in_cgs = r_200_cgs
-
-print("UnitMass_in_cgs: ", const_unit_mass_in_cgs)
-print("UnitLength_in_cgs: ", const_unit_length_in_cgs)
-print("UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs)
-
-# Now set the magnetic field unit
-
-const_unit_magnetic_field_in_cgs = 1e-7 # 1muG
-
-# Derived quantities
-const_unit_time_in_cgs = const_unit_length_in_cgs / const_unit_velocity_in_cgs
-const_unit_current_in_cgs = const_unit_mass_in_cgs / (
- const_unit_magnetic_field_in_cgs * const_unit_time_in_cgs ** 2
-)
-
-print("UnitTime_in_cgs: ", const_unit_time_in_cgs)
-print("UnitCurrent_in_cgs: ", const_unit_current_in_cgs)
-
-const_G = (
- CONST_G_CGS
- * const_unit_mass_in_cgs
- * const_unit_time_in_cgs
- * const_unit_time_in_cgs
- / (const_unit_length_in_cgs * const_unit_length_in_cgs * const_unit_length_in_cgs)
-)
-print("G=", const_G)
-
-# Parameters
-periodic = 1 # 1 For periodic box
-boxSize = 4.0
-G = const_G
-N = int(sys.argv[1]) # Number of particles
-
-# Create the file
-filename = "CoolingHalo.hdf5"
-file = h5py.File(filename, "w")
-
-# Units
-grp = file.create_group("/Units")
-grp.attrs["Unit length in cgs (U_L)"] = const_unit_length_in_cgs
-grp.attrs["Unit mass in cgs (U_M)"] = const_unit_mass_in_cgs
-grp.attrs["Unit time in cgs (U_t)"] = (
- const_unit_length_in_cgs / const_unit_velocity_in_cgs
-)
-grp.attrs["Unit current in cgs (U_I)"] = const_unit_current_in_cgs
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
-
-# set seed for random number
-np.random.seed(1234)
-
-gas_mass = (
- f_b * np.sqrt(3.0) / 2.0
-) # virial mass of halo is 1, virial radius is 1, enclosed mass scales with r
-gas_particle_mass = gas_mass / float(N)
-
-# Positions
-# r^(-2) distribution corresponds to uniform distribution in radius
-radius = (
- boxSize * np.sqrt(3.0) / 2.0 * np.random.rand(N)
-) # the diagonal extent of the cube
-ctheta = -1.0 + 2 * np.random.rand(N)
-stheta = np.sqrt(1.0 - ctheta ** 2)
-phi = 2 * math.pi * np.random.rand(N)
-coords = np.zeros((N, 3))
-coords[:, 0] = radius * stheta * np.cos(phi)
-coords[:, 1] = radius * stheta * np.sin(phi)
-coords[:, 2] = radius * ctheta
-
-# shift to centre of box
-coords += np.full((N, 3), boxSize / 2.0)
-print("x range = (%f,%f)" % (np.min(coords[:, 0]), np.max(coords[:, 0])))
-print("y range = (%f,%f)" % (np.min(coords[:, 1]), np.max(coords[:, 1])))
-print("z range = (%f,%f)" % (np.min(coords[:, 2]), np.max(coords[:, 2])))
-
-print(np.mean(coords[:, 0]))
-print(np.mean(coords[:, 1]))
-print(np.mean(coords[:, 2]))
-
-# now find the particles which are within the box
-
-x_coords = coords[:, 0]
-y_coords = coords[:, 1]
-z_coords = coords[:, 2]
-
-ind = np.where(x_coords < boxSize)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-ind = np.where(x_coords > 0.0)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-ind = np.where(y_coords < boxSize)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-ind = np.where(y_coords > 0.0)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-ind = np.where(z_coords < boxSize)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-ind = np.where(z_coords > 0.0)[0]
-x_coords = x_coords[ind]
-y_coords = y_coords[ind]
-z_coords = z_coords[ind]
-
-# count number of particles
-
-N = x_coords.size
-
-print("Number of particles in the box = ", N)
-
-# make the coords and radius arrays again
-coords = np.zeros((N, 3))
-coords[:, 0] = x_coords
-coords[:, 1] = y_coords
-coords[:, 2] = z_coords
-
-radius = np.sqrt(
- (coords[:, 0] - boxSize / 2.0) ** 2
- + (coords[:, 1] - boxSize / 2.0) ** 2
- + (coords[:, 2] - boxSize / 2.0) ** 2
-)
-
-# now give particle's velocities and magnetic fields
-v = np.zeros((N, 3))
-B = np.zeros((N, 3))
-
-# first work out total angular momentum of the halo within the virial radius
-# we work in units where r_vir = 1 and M_vir = 1
-Total_E = v_200 ** 2 / 2.0
-J = spin_lambda * const_G / np.sqrt(Total_E)
-print("J =", J)
-# all particles within the virial radius have omega parallel to the z-axis, magnitude
-# is proportional to 1 over the radius
-omega = np.zeros((N, 3))
-for i in range(N):
- omega[i, 2] = 3.0 * J / radius[i]
- v[i, :] = np.cross(omega[i, :], (coords[i, :] - boxSize / 2.0))
- B[i, 2] = B0_cgs / const_unit_magnetic_field_in_cgs
- #B[i, 0] = B0_cgs / const_unit_magnetic_field_in_cgs
-
-# Header
-grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = boxSize
-grp.attrs["NumPart_Total"] = [N, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_ThisFile"] = [N, 0, 0, 0, 0, 0]
-grp.attrs["Time"] = 0.0
-grp.attrs["NumFilesPerSnapshot"] = 1
-grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["Dimension"] = 3
-
-# Particle group
-grp = file.create_group("/PartType0")
-
-ds = grp.create_dataset("Coordinates", (N, 3), "d")
-ds[()] = coords
-coords = np.zeros(1)
-
-ds = grp.create_dataset("Velocities", (N, 3), "f")
-ds[()] = v
-v = np.zeros(1)
-
-# Magnetic field
-ds = grp.create_dataset("MagneticFluxDensities", (N, 3), "f")
-ds[()] = B
-B = np.zeros(1)
-
-# All particles of equal mass
-m = np.full((N,), gas_particle_mass)
-ds = grp.create_dataset("Masses", (N,), "f")
-ds[()] = m
-m = np.zeros(1)
-
-# Smoothing lengths
-l = (4.0 * np.pi * radius ** 2 / N) ** (
- 1.0 / 3.0
-) # local mean inter-particle separation
-h = np.full((N,), eta * l)
-ds = grp.create_dataset("SmoothingLength", (N,), "f")
-ds[()] = h
-h = np.zeros(1)
-
-# Internal energies
-u = v_200 ** 2 / (2.0 * (gamma - 1.0))
-u = np.full((N,), u)
-ds = grp.create_dataset("InternalEnergy", (N,), "f")
-ds[()] = u
-u = np.zeros(1)
-
-# Particle IDs
-ids = 1 + np.linspace(0, N, N, endpoint=False)
-ds = grp.create_dataset("ParticleIDs", (N,), "L")
-ds[()] = ids
-
-file.close()
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/particle_tracker.py b/examples/MHDTests/CoolingHaloWithSpin_MHD/particle_tracker.py
index 4da3967ab95f54cb1bee61120d3c784304b16092..20ed3487ee4a41bfe6aafa75b869be18f5ea8583 100644
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/particle_tracker.py
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/particle_tracker.py
@@ -6,43 +6,59 @@ import matplotlib.pyplot as plt
from matplotlib import rc
+
# Set the global font to be DejaVu Sans, size 10 (or any other sans-seri f font of your choice!)
-rc('font',**{'family':'sans-serif','sans-serif':['DejaVu Sans'],'size':10})
+rc("font", **{"family": "sans-serif", "sans-serif": ["DejaVu Sans"], "size": 10})
# Set the font used for MathJax - more on this later
-rc('mathtext',**{'default':'regular'})
-
-
-
-# setting units here
-units_dict = {'length':3.08e21*u.cm,'density':1.98e33*u.g/(3.08e21*u.cm)**3,'velocity': u.km/u.s,'momentum': 1.98e33*u.g*u.km/u.s, 'pressure': 1e10*(1.98e33*u.g)/((3.08e21*u.cm)*(3.15e16 * u.s)**2), 'magneticfield':1e-7*u.g / (u.s**2 * u.statA),'magneticfieldpertime':1e-7*u.g/(u.s**2 * u.statA * 3.156e7 * u.s),'magneticfieldgradient':1e-7*u.g/(u.s**2 * u.statA * u.cm), 'time':3.15e16 * u.s, 'mass':1.98e33*u.g,'acceleration':u.cm/u.s**2}
-#parsec - 3.086e18*u.cm
+rc("mathtext", **{"default": "regular"})
+
+
+# setting units here
+units_dict = {
+ "length": 3.08e21 * u.cm,
+ "density": 1.98e33 * u.g / (3.08e21 * u.cm) ** 3,
+ "velocity": u.km / u.s,
+ "momentum": 1.98e33 * u.g * u.km / u.s,
+ "pressure": 1e10 * (1.98e33 * u.g) / ((3.08e21 * u.cm) * (3.15e16 * u.s) ** 2),
+ "magneticfield": 1e-7 * u.g / (u.s ** 2 * u.statA),
+ "magneticfieldpertime": 1e-7 * u.g / (u.s ** 2 * u.statA * 3.156e7 * u.s),
+ "magneticfieldgradient": 1e-7 * u.g / (u.s ** 2 * u.statA * u.cm),
+ "time": 3.15e16 * u.s,
+ "mass": 1.98e33 * u.g,
+ "acceleration": u.cm / u.s ** 2,
+}
+# parsec - 3.086e18*u.cm
# free fall time estimate
# Constants
-G = 6.67430e-11 * u.m **3 / u.kg / u.s **2 #6.67430e-8
-G = G.to(u.cm **3 / u.g / u.s **2)
+G = 6.67430e-11 * u.m ** 3 / u.kg / u.s ** 2 # 6.67430e-8
+G = G.to(u.cm ** 3 / u.g / u.s ** 2)
# Parameters (taken from Hopkins 2016)
Rcloud = 4.628516371e16 * u.cm
M = 1.99e33 * u.g # total mass of the sphere
-rhocloud0 = M/(4/3 * np.pi * Rcloud**3 )
-rhouniform = rhocloud0/360
-t_ff_lit = np.sqrt(3/(2*np.pi*G*rhocloud0))
-t_ff_wiki = np.sqrt(3*np.pi/(32*G*rhocloud0))
+rhocloud0 = M / (4 / 3 * np.pi * Rcloud ** 3)
+rhouniform = rhocloud0 / 360
+t_ff_lit = np.sqrt(3 / (2 * np.pi * G * rhocloud0))
+t_ff_wiki = np.sqrt(3 * np.pi / (32 * G * rhocloud0))
t_ff_script = 3e4 * 3.156e7 * u.s
+
def abs_vec(vec):
res = np.sqrt(vec[:, 0] ** 2 + vec[:, 1] ** 2 + vec[:, 2] ** 2)
return res
+
def dot_vec(vec1, vec2):
- res = vec1[:, 0] * vec2[:, 0] + vec1[:, 1] * vec2[:, 1] + vec1[:, 2]* vec2[:, 2]
+ res = vec1[:, 0] * vec2[:, 0] + vec1[:, 1] * vec2[:, 1] + vec1[:, 2] * vec2[:, 2]
return res
-def get_projection(vec,direction_vec):
- res = dot_vec(vec,direction_vec)/abs_vec(direction_vec)
+
+def get_projection(vec, direction_vec):
+ res = dot_vec(vec, direction_vec) / abs_vec(direction_vec)
return res
+
def cross_vec(vec1, vec2):
res_vec = np.zeros((len(vec1), 3))
res_vec[:, 0] = vec1[:, 1] * vec2[:, 2] - vec1[:, 2] * vec2[:, 1]
@@ -54,7 +70,7 @@ def cross_vec(vec1, vec2):
def extract_variables_of_interest(data):
# set constant
- mu0 = 1.25663706127e-1 * u.g * u.cm / (u.s**2 * u.statA**2)
+ mu0 = 1.25663706127e-1 * u.g * u.cm / (u.s ** 2 * u.statA ** 2)
# Get snapshot parameters
time = data.metadata.time
@@ -64,57 +80,88 @@ def extract_variables_of_interest(data):
pids = data.gas.particle_ids
# Get physical quantities
- coordinates = data.gas.coordinates.to_physical().astype('float64')
- densities = data.gas.densities.to_physical().astype('float64')
- smoothing_lengths = data.gas.smoothing_lengths.to_physical().astype('float64')
- velocities = data.gas.velocities.to_physical().astype('float64')
- pressures = data.gas.pressures.to_physical().astype('float64')
- masses = data.gas.masses.to_physical().astype('float64')
+ coordinates = data.gas.coordinates.to_physical().astype("float64")
+ densities = data.gas.densities.to_physical().astype("float64")
+ smoothing_lengths = data.gas.smoothing_lengths.to_physical().astype("float64")
+ velocities = data.gas.velocities.to_physical().astype("float64")
+ pressures = data.gas.pressures.to_physical().astype("float64")
+ masses = data.gas.masses.to_physical().astype("float64")
momentums = velocities * masses[:, np.newaxis]
- magnetic_flux_densities = data.gas.magnetic_flux_densities.to_physical().astype('float64')
- magnetic_flux_densitiesdt = data.gas.magnetic_flux_densitiesdt.to_physical().astype('float64')
- magnetic_divergences = data.gas.magnetic_divergences.to_physical().astype('float64')
+ magnetic_flux_densities = data.gas.magnetic_flux_densities.to_physical().astype(
+ "float64"
+ )
+ magnetic_flux_densitiesdt = data.gas.magnetic_flux_densitiesdt.to_physical().astype(
+ "float64"
+ )
+ magnetic_divergences = data.gas.magnetic_divergences.to_physical().astype("float64")
- magnetic_pressures = dot_vec(magnetic_flux_densities,magnetic_flux_densities)/(2*mu0)
- dynamic_pressures = densities*dot_vec(velocities, velocities)
+ magnetic_pressures = dot_vec(magnetic_flux_densities, magnetic_flux_densities) / (
+ 2 * mu0
+ )
+ dynamic_pressures = densities * dot_vec(velocities, velocities)
# Center of mass position tracking
CM_frame_coordinates = coordinates
for k in range(3):
- CM_k = np.sum(coordinates[:,k].value * masses.value)/np.sum(masses.value)*coordinates.units
+ CM_k = (
+ np.sum(coordinates[:, k].value * masses.value)
+ / np.sum(masses.value)
+ * coordinates.units
+ )
CM_frame_coordinates -= CM_k
# working with units:
- coordinates = coordinates.to(units_dict['length'])
- CM_frame_coordinates = CM_frame_coordinates.to(units_dict['length'])
- densities = densities.to(units_dict['density'])
- smoothing_lengths = smoothing_lengths.to(units_dict['length'])
- velocities = velocities.to(units_dict['velocity'])
- momentums = momentums.to(units_dict['momentum'])
- pressures = pressures.to(units_dict['pressure'])
- magnetic_flux_densitiesdt = magnetic_flux_densitiesdt.to(units_dict['magneticfieldpertime'])
- magnetic_divergences = magnetic_divergences.to(units_dict['magneticfieldgradient'])
- masses = masses.to(units_dict['mass'])
- magnetic_pressures = magnetic_pressures.to(units_dict['pressure'])
- dynamic_pressures = dynamic_pressures.to(units_dict['pressure'])
+ coordinates = coordinates.to(units_dict["length"])
+ CM_frame_coordinates = CM_frame_coordinates.to(units_dict["length"])
+ densities = densities.to(units_dict["density"])
+ smoothing_lengths = smoothing_lengths.to(units_dict["length"])
+ velocities = velocities.to(units_dict["velocity"])
+ momentums = momentums.to(units_dict["momentum"])
+ pressures = pressures.to(units_dict["pressure"])
+ magnetic_flux_densitiesdt = magnetic_flux_densitiesdt.to(
+ units_dict["magneticfieldpertime"]
+ )
+ magnetic_divergences = magnetic_divergences.to(units_dict["magneticfieldgradient"])
+ masses = masses.to(units_dict["mass"])
+ magnetic_pressures = magnetic_pressures.to(units_dict["pressure"])
+ dynamic_pressures = dynamic_pressures.to(units_dict["pressure"])
magnetic_flux_densities_norm = abs_vec(magnetic_flux_densities)
- time = time.to(units_dict['time'])
+ time = time.to(units_dict["time"])
+
+ r0_no_cut = (
+ magnetic_divergences
+ * smoothing_lengths
+ / (abs_vec(magnetic_flux_densities.value) * magnetic_flux_densities.units)
+ )
- r0_no_cut = magnetic_divergences * smoothing_lengths / (abs_vec(magnetic_flux_densities.value)*magnetic_flux_densities.units)
+ CM_frame_z = np.abs(CM_frame_coordinates[:, 2])
+ data_dict = {
+ "pids": pids,
+ "coordinates": coordinates,
+ "masses": masses,
+ "densities": densities,
+ "smoothing_lengths": smoothing_lengths,
+ "velocities": velocities,
+ "momentums": momentums,
+ "pressures": pressures,
+ "magnetic_flux_densities": magnetic_flux_densities,
+ "CM_frame_coordinates": CM_frame_coordinates,
+ "CM_frame_z": CM_frame_z,
+ "magnetic_flux_densities_norm": magnetic_flux_densities_norm,
+ "magnetic_pressures": magnetic_pressures,
+ "dynamic_pressures": dynamic_pressures,
+ }
- CM_frame_z = np.abs(CM_frame_coordinates[:,2])
- data_dict = {'pids':pids,'coordinates':coordinates,'masses':masses,'densities':densities,'smoothing_lengths':smoothing_lengths,'velocities':velocities,'momentums':momentums,'pressures':pressures,'magnetic_flux_densities':magnetic_flux_densities,'CM_frame_coordinates':CM_frame_coordinates,'CM_frame_z':CM_frame_z,'magnetic_flux_densities_norm':magnetic_flux_densities_norm,'magnetic_pressures':magnetic_pressures,'dynamic_pressures':dynamic_pressures}
-
values = {}
units = {}
for key in data_dict.keys():
- values = values | {key:data_dict[key].value}
+ values = values | {key: data_dict[key].value}
- snapshot_parameters = {'time':time.value, 'z':z, 'a':a}
+ snapshot_parameters = {"time": time.value, "z": z, "a": a}
# Retrieve some information about the simulation run
artDiffusion = data.metadata.hydro_scheme["Artificial Diffusion Constant"]
@@ -128,222 +175,277 @@ def extract_variables_of_interest(data):
kernel = data.metadata.hydro_scheme["Kernel function"]
neighbours = data.metadata.hydro_scheme["Kernel target N_ngb"]
- metadata = {'artDiffusion':artDiffusion, 'dedHyp':dedHyp, 'dedHypDivv':dedHypDivv, 'dedPar':dedPar,'eta':eta,'git':git,'gitBranch':gitBranch,'hydroScheme':hydroScheme,'kernel':kernel,'neighbours':neighbours}
+ metadata = {
+ "artDiffusion": artDiffusion,
+ "dedHyp": dedHyp,
+ "dedHypDivv": dedHypDivv,
+ "dedPar": dedPar,
+ "eta": eta,
+ "git": git,
+ "gitBranch": gitBranch,
+ "hydroScheme": hydroScheme,
+ "kernel": kernel,
+ "neighbours": neighbours,
+ }
+
+ return {"values": values, "parameters": snapshot_parameters, "metadata": metadata}
- return {'values':values,'parameters':snapshot_parameters,'metadata':metadata}
def updoad_shapshots_data_from_folder(folderaddr):
- addr_book = glob(folderaddr+'/**/*.hdf5',recursive=True)
+ addr_book = glob(folderaddr + "/**/*.hdf5", recursive=True)
addr_book = sorted(addr_book)
- names = [addr.split('/')[-1].split('.hdf5')[0] for addr in addr_book]
- #print(names)
+ names = [addr.split("/")[-1].split(".hdf5")[0] for addr in addr_book]
+ # print(names)
snapshots_data = []
for addr in addr_book:
snapshot = load(addr)
- snapshots_data+=[extract_variables_of_interest(snapshot)]
+ snapshots_data += [extract_variables_of_interest(snapshot)]
return snapshots_data, names
-def upload_particle_history(all_data,pid):
+
+def upload_particle_history(all_data, pid):
particle_data = []
for snapshot in all_data:
- particle_snapshot_dict = snapshot['parameters']
- indx = np.argwhere(snapshot['values']['pids']==pid)[0][0]
- for key in snapshot['values'].keys():
- particle_snapshot_dict = particle_snapshot_dict | {key:snapshot['values'][key][indx]}
- particle_data+=[particle_snapshot_dict]
-
- # transform list of dicts to single
+ particle_snapshot_dict = snapshot["parameters"]
+ indx = np.argwhere(snapshot["values"]["pids"] == pid)[0][0]
+ for key in snapshot["values"].keys():
+ particle_snapshot_dict = particle_snapshot_dict | {
+ key: snapshot["values"][key][indx]
+ }
+ particle_data += [particle_snapshot_dict]
+
+ # transform list of dicts to single
# Dynamically extract keys from the first dictionary
keys = particle_data[0].keys()
# Transform the list of dictionaries to a single dictionary
result = {key: np.array([entry[key] for entry in particle_data]) for key in keys}
- result = result | {'pid':pid}
+ result = result | {"pid": pid}
return result
-def identify_largest_quantity_particles(all_snapshots,quantity,isvec=False,region_cut=True):
+
+def identify_largest_quantity_particles(
+ all_snapshots, quantity, isvec=False, region_cut=True
+):
largest_quantity_particles = []
for snap in all_snapshots:
- quantity_values = snap['values'][quantity]
- zcoords = np.abs(snap['values']['CM_frame_coordinates'][:,2])
- rcoords = np.sqrt(snap['values']['CM_frame_coordinates'][:,0]**2+snap['values']['CM_frame_coordinates'][:,1]**2)
- velocities = snap['values']['velocities']
- vabs = np.sqrt(velocities[:,0]**2+velocities[:,1]**2+velocities[:,2]**2)
+ quantity_values = snap["values"][quantity]
+ zcoords = np.abs(snap["values"]["CM_frame_coordinates"][:, 2])
+ rcoords = np.sqrt(
+ snap["values"]["CM_frame_coordinates"][:, 0] ** 2
+ + snap["values"]["CM_frame_coordinates"][:, 1] ** 2
+ )
+ velocities = snap["values"]["velocities"]
+ vabs = np.sqrt(
+ velocities[:, 0] ** 2 + velocities[:, 1] ** 2 + velocities[:, 2] ** 2
+ )
mask = np.array([True for i in range(len(quantity_values))])
if region_cut:
- zcut = (zcoords > 0.015*0.015) & (zcoords < 0.15*0.015) # z>0.015 Rcloud
- rcut = (rcoords < 0.15*0.015)
-
+ zcut = (zcoords > 0.015 * 0.015) & (
+ zcoords < 0.15 * 0.015
+ ) # z>0.015 Rcloud
+ rcut = rcoords < 0.15 * 0.015
+
mask = (zcut) & (rcut)
- #quantity_values = quantity_values[zcut]
+ # quantity_values = quantity_values[zcut]
if isvec:
- quantity_values = np.sqrt(quantity_values[:,0]**2+quantity_values[:,1]**2+quantity_values[:,2]**2)
- pids = snap['values']['pids']
+ quantity_values = np.sqrt(
+ quantity_values[:, 0] ** 2
+ + quantity_values[:, 1] ** 2
+ + quantity_values[:, 2] ** 2
+ )
+ pids = snap["values"]["pids"]
quantity_values_cut = quantity_values[mask]
- #print('snapshot at t=',snap['parameters']['time'])
+ # print('snapshot at t=',snap['parameters']['time'])
indx = np.argmax(quantity_values_cut)
- #print('particle #',pids[mask][indx])
- #print('quantity value = ',quantity_values[pids[mask][indx]])
+ # print('particle #',pids[mask][indx])
+ # print('quantity value = ',quantity_values[pids[mask][indx]])
largest_quantity_particles.append(pids[mask][indx])
return largest_quantity_particles
-def plot_quatities_for_particle_vs_time(particle_history,quantity_list,outputfileaddr='./quantities.png',time_var = 'time',timeinterval=[0,3],title='quantities vs time',customtimeinterval = False,customquantityinterval = False, quantityinterval=[1e-1,1e1],logscale=True):
- fig, ax = plt.subplots(1, len(quantity_list), figsize=(6.5*len(quantity_list), 6))
+def plot_quatities_for_particle_vs_time(
+ particle_history,
+ quantity_list,
+ outputfileaddr="./quantities.png",
+ time_var="time",
+ timeinterval=[0, 3],
+ title="quantities vs time",
+ customtimeinterval=False,
+ customquantityinterval=False,
+ quantityinterval=[1e-1, 1e1],
+ logscale=True,
+):
+
+ fig, ax = plt.subplots(1, len(quantity_list), figsize=(6.5 * len(quantity_list), 6))
for i in range(len(quantity_list)):
quantity = particle_history[quantity_list[i]]
times = particle_history[time_var]
- if len(quantity.shape)>1:
- quantity = np.sqrt(quantity[:,0]**2+quantity[:,1]**2+quantity[:,2]**2)
-
- mask_nonzero = quantity!=0.0
- #quantity=quantity[mask_nonzero]
- #times = times[mask_nonzero]
- mask_in_range = (times>timeinterval[0]) & (times<timeinterval[-1])
- #print(len(times[mask_in_range]))
- #if len(quantity[mask_in_range])!=0:
- #ax[i].axvline(x=1.378579,ymin = 0.1*np.min(quantity[mask_in_range]),ymax = 10*np.max(quantity[mask_in_range]),color='red',linestyle='dashed')
- #ax[i].axvline(x=1.378579,ymin = 1e-30,ymax = 1e30,color='red',linestyle='dashed')
-
- #handle only positive and positive-negative values:
+ if len(quantity.shape) > 1:
+ quantity = np.sqrt(
+ quantity[:, 0] ** 2 + quantity[:, 1] ** 2 + quantity[:, 2] ** 2
+ )
+
+ mask_nonzero = quantity != 0.0
+ # quantity=quantity[mask_nonzero]
+ # times = times[mask_nonzero]
+ mask_in_range = (times > timeinterval[0]) & (times < timeinterval[-1])
+ # print(len(times[mask_in_range]))
+ # if len(quantity[mask_in_range])!=0:
+ # ax[i].axvline(x=1.378579,ymin = 0.1*np.min(quantity[mask_in_range]),ymax = 10*np.max(quantity[mask_in_range]),color='red',linestyle='dashed')
+ # ax[i].axvline(x=1.378579,ymin = 1e-30,ymax = 1e30,color='red',linestyle='dashed')
+
+ # handle only positive and positive-negative values:
if logscale:
- mask = quantity<0
- if np.sum(mask)>0:
+ mask = quantity < 0
+ if np.sum(mask) > 0:
# Separate positive and negative values
- positive = np.where(quantity > 0, quantity, np.nan) # Use NaN for negative values
- negative = np.where(quantity < 0, quantity, np.nan) # Use NaN for positive values
- ax[i].plot(times, positive,color='red',marker='.')
- ax[i].plot(times,-negative,color='blue',marker='.')
+ positive = np.where(
+ quantity > 0, quantity, np.nan
+ ) # Use NaN for negative values
+ negative = np.where(
+ quantity < 0, quantity, np.nan
+ ) # Use NaN for positive values
+ ax[i].plot(times, positive, color="red", marker=".")
+ ax[i].plot(times, -negative, color="blue", marker=".")
else:
- ax[i].plot(times,quantity,color='black',marker='.')
+ ax[i].plot(times, quantity, color="black", marker=".")
else:
- ax[i].plot(times,quantity,color='black',marker='.')
- ax[i].set_ylabel(quantity_list[i],fontsize=20)
- ax[i].set_xlabel('$t/t_{ff}$',fontsize=20)
+ ax[i].plot(times, quantity, color="black", marker=".")
+ ax[i].set_ylabel(quantity_list[i], fontsize=20)
+ ax[i].set_xlabel("$t/t_{ff}$", fontsize=20)
if logscale:
- ax[i].set_yscale('log')
- ax[i].tick_params(axis='both', labelsize=20)
+ ax[i].set_yscale("log")
+ ax[i].tick_params(axis="both", labelsize=20)
if customtimeinterval:
ax[i].set_xlim(timeinterval)
if customquantityinterval:
ax[i].set_ylim(quantityinterval)
ax[i].grid()
- fig.suptitle('Particle #'+str(particle_history['pid']),fontsize=20)
+ fig.suptitle("Particle #" + str(particle_history["pid"]), fontsize=20)
fig.tight_layout()
plt.savefig(outputfileaddr)
plt.close(fig)
-
return 0
-Lcut = 15 # Lcut in kPc.
-def plot_quantity_vs_time(all_snapshots,quantity,outputfileaddr,isvec=False,region_cut=True,label=''):
+
+Lcut = 15 # Lcut in kPc.
+
+
+def plot_quantity_vs_time(
+ all_snapshots, quantity, outputfileaddr, isvec=False, region_cut=True, label=""
+):
quantity_array = []
rmsq_array = []
mzq_array = []
times_array = []
npart_array = []
for snap in all_snapshots:
- quantity_values = snap['values'][quantity]
+ quantity_values = snap["values"][quantity]
Np = len(quantity_values)
- #total_quantity = np.sum(quantity_values,axis=0)
- #rms_quantity = np.sqrt(np.mean(quantity_values[:,0]**2+quantity_values[:,1]**2+quantity_values[:,2]**2))
- print(len(quantity_values))
+ # total_quantity = np.sum(quantity_values,axis=0)
+ # rms_quantity = np.sqrt(np.mean(quantity_values[:,0]**2+quantity_values[:,1]**2+quantity_values[:,2]**2))
+ print(len(quantity_values))
if region_cut:
- z = snap['values']['CM_frame_coordinates'][:,2]
- y = snap['values']['CM_frame_coordinates'][:,1]
- x = snap['values']['CM_frame_coordinates'][:,0]
- mask = ((np.abs(z)<=Lcut) & (np.abs(y)<=Lcut) & (np.abs(x)<=Lcut))
+ z = snap["values"]["CM_frame_coordinates"][:, 2]
+ y = snap["values"]["CM_frame_coordinates"][:, 1]
+ x = snap["values"]["CM_frame_coordinates"][:, 0]
+ mask = (np.abs(z) <= Lcut) & (np.abs(y) <= Lcut) & (np.abs(x) <= Lcut)
quantity_values = quantity_values[mask]
- total_quantity = np.sum(quantity_values,axis=0)
- rms_quantity = np.sqrt(np.mean(quantity_values[:,0]**2+quantity_values[:,1]**2+quantity_values[:,2]**2))
- mz_quantity = np.mean(np.abs(quantity_values[:,2]))
- time = snap['parameters']['time']
+ total_quantity = np.sum(quantity_values, axis=0)
+ rms_quantity = np.sqrt(
+ np.mean(
+ quantity_values[:, 0] ** 2
+ + quantity_values[:, 1] ** 2
+ + quantity_values[:, 2] ** 2
+ )
+ )
+ mz_quantity = np.mean(np.abs(quantity_values[:, 2]))
+ time = snap["parameters"]["time"]
quantity_array.append(total_quantity)
rmsq_array.append(rms_quantity)
mzq_array.append(mz_quantity)
- times_array.append(time)
+ times_array.append(time)
npart_array.append(len(quantity_values))
-
# Retrieve some information about the simulation run
- artDiffusion = snap['metadata']['artDiffusion']
- dedHyp = snap['metadata']['dedHyp']
- dedHypDivv = snap['metadata']['dedHypDivv']
- dedPar = snap['metadata']['dedPar']
- eta = snap['metadata']['eta']
- git = snap['metadata']['git']
- gitBranch = snap['metadata']['gitBranch']
- hydroScheme = snap['metadata']['hydroScheme']
- kernel = snap['metadata']['kernel']
- neighbours = snap['metadata']['neighbours']
-
-
+ artDiffusion = snap["metadata"]["artDiffusion"]
+ dedHyp = snap["metadata"]["dedHyp"]
+ dedHypDivv = snap["metadata"]["dedHypDivv"]
+ dedPar = snap["metadata"]["dedPar"]
+ eta = snap["metadata"]["eta"]
+ git = snap["metadata"]["git"]
+ gitBranch = snap["metadata"]["gitBranch"]
+ hydroScheme = snap["metadata"]["hydroScheme"]
+ kernel = snap["metadata"]["kernel"]
+ neighbours = snap["metadata"]["neighbours"]
+
quantity_array = np.array(quantity_array)
rmsq_array = np.array(rmsq_array)
mzq_array = np.array(mzq_array)
times_array = np.array(times_array)
npart_array = np.array(npart_array)
-
- fig, ax = plt.subplots(1, 3, figsize=(5.5*3, 5))
+ fig, ax = plt.subplots(1, 3, figsize=(5.5 * 3, 5))
if isvec:
- quantity_abs = np.sqrt(quantity_array[:,0]**2+quantity_array[:,1]**2+quantity_array[:,2]**2)
+ quantity_abs = np.sqrt(
+ quantity_array[:, 0] ** 2
+ + quantity_array[:, 1] ** 2
+ + quantity_array[:, 2] ** 2
+ )
else:
quantity_abs = np.abs(quantity_array)
- #ax[0].plot(times_array,quantity_array[:,2],color='black',marker='.', label='$p_z$')
- #ax[0].plot(times_array,quantity_array[:,0],color='blue' ,marker='.', label='$p_x$')
- #ax[0].plot(times_array,quantity_array[:,1],color='red' ,marker='.', label='$p_y$')
- #ax[0].plot(times_array,np.abs(quantity_array[:,2]),color='black',marker='.', label='$p_z$')
- #pq = quantity_array[:,2]#/npart_array
- #mask = pq<0
- #lpq = '$mean p_z$'
- #if np.sum(mask)>0:
+ # ax[0].plot(times_array,quantity_array[:,2],color='black',marker='.', label='$p_z$')
+ # ax[0].plot(times_array,quantity_array[:,0],color='blue' ,marker='.', label='$p_x$')
+ # ax[0].plot(times_array,quantity_array[:,1],color='red' ,marker='.', label='$p_y$')
+ # ax[0].plot(times_array,np.abs(quantity_array[:,2]),color='black',marker='.', label='$p_z$')
+ # pq = quantity_array[:,2]#/npart_array
+ # mask = pq<0
+ # lpq = '$mean p_z$'
+ # if np.sum(mask)>0:
# # Separate positive and negative values
# positive = np.where(pq > 0, pq, np.nan) # Use NaN for negative values
# negative = np.where(pq < 0, pq, np.nan) # Use NaN for positive values
# ax[0].plot(times_array, positive,color='red', marker='.',label = 'mean $p_z>0$')
# ax[0].plot(times_array,-negative,color='blue',marker='.',label = 'mean $p_z<0$')
- #else:
- #ax[0].plot(times_array,quantity_abs,color='black',marker='.',label = 'B')
+ # else:
+ # ax[0].plot(times_array,quantity_abs,color='black',marker='.',label = 'B')
- ax[0].plot(times_array, rmsq_array, color='green', marker='.',label='$B_{rms}$')
- #ax[0].plot(times_array, mzq_array, color='black', marker='.',label='mean $|p_z|$')
+ ax[0].plot(times_array, rmsq_array, color="green", marker=".", label="$B_{rms}$")
+ # ax[0].plot(times_array, mzq_array, color='black', marker='.',label='mean $|p_z|$')
+
+ # ax[0].plot(times_array,np.abs(quantity_array[:,0]),color='blue' ,marker='.', label='$p_x$')
+ # ax[0].plot(times_array,np.abs(quantity_array[:,1]),color='red' ,marker='.', label='$p_y$')
- #ax[0].plot(times_array,np.abs(quantity_array[:,0]),color='blue' ,marker='.', label='$p_x$')
- #ax[0].plot(times_array,np.abs(quantity_array[:,1]),color='red' ,marker='.', label='$p_y$')
-
ax[0].legend()
ax[0].grid()
- #ax[0].set_yscale('log')
- ax[0].set_ylabel('$|B|$, [$\mu G$]',fontsize=20)
- ax[0].set_xlabel('$time [Gyr]$',fontsize=20)
- ax[0].set_ylim([0,10.0])
- ax[0].set_xlim([0,3.0])
- #ax[0].set_ylim([1e-8,1e-4])
-
- ax[1].plot(times_array,npart_array,color='black',marker='.')
+ # ax[0].set_yscale('log')
+ ax[0].set_ylabel("$|B|$, [$\mu G$]", fontsize=20)
+ ax[0].set_xlabel("$time [Gyr]$", fontsize=20)
+ ax[0].set_ylim([0, 10.0])
+ ax[0].set_xlim([0, 3.0])
+ # ax[0].set_ylim([1e-8,1e-4])
+
+ ax[1].plot(times_array, npart_array, color="black", marker=".")
ax[1].grid()
- ax[1].set_yscale('log')
- ax[1].set_ylabel(f'$N_p$'+label ,fontsize=20) #, Lcut = {Lcut} A.U.',fontsize=20)
- ax[1].set_xlabel('$time [Gyr]$',fontsize=20)
- ax[1].set_ylim([1e1,1e5])
- ax[1].set_xlim([0,3.0])
+ ax[1].set_yscale("log")
+ ax[1].set_ylabel(
+ f"$N_p$" + label, fontsize=20
+ ) # , Lcut = {Lcut} A.U.',fontsize=20)
+ ax[1].set_xlabel("$time [Gyr]$", fontsize=20)
+ ax[1].set_ylim([1e1, 1e5])
+ ax[1].set_xlim([0, 3.0])
# add panel with infromation about the run
text_common_args = dict(
fontsize=10, ha="center", va="center", transform=ax[2].transAxes
)
- ax[2].text(
- 0.5,
- 0.8,
- "Cooling halo with spin, $Np=%.0f$ " % Np,
- **text_common_args,
- )
+ ax[2].text(0.5, 0.8, "Cooling halo with spin, $Np=%.0f$ " % Np, **text_common_args)
ax[2].text(0.5, 0.7, "swift %s" % git.decode("utf-8"), **text_common_args)
ax[2].text(0.5, 0.6, "Branch %s" % gitBranch.decode("utf-8"), **text_common_args)
ax[2].text(0.5, 0.5, hydroScheme.decode("utf-8"), **text_common_args)
@@ -367,36 +469,47 @@ def plot_quantity_vs_time(all_snapshots,quantity,outputfileaddr,isvec=False,regi
)
ax[2].axis("off")
-
-
- #fig.suptitle(quantity,fontsize=20)
+ # fig.suptitle(quantity,fontsize=20)
fig.tight_layout()
plt.savefig(outputfileaddr)
plt.close(fig)
- return
-
-def plot_pressure_vs_time(all_snapshots,outputfileaddr,isvec=False,region_cut=True,label='',rcut=15,zcut=0.5):
+ return
+
+
+def plot_pressure_vs_time(
+ all_snapshots,
+ outputfileaddr,
+ isvec=False,
+ region_cut=True,
+ label="",
+ rcut=15,
+ zcut=0.5,
+):
pressures_array = []
magnetic_pressures_array = []
dynamic_pressures_array = []
times_array = []
npart_array = []
- rcut*=units_dict['length']
- zcut*=units_dict['length']
+ rcut *= units_dict["length"]
+ zcut *= units_dict["length"]
for snap in all_snapshots:
- pressures_values = snap['values']['pressures']*units_dict['pressure']
- magnetic_pressures_values = snap['values']['magnetic_pressures']*units_dict['pressure']
- dynamic_pressures_values = snap['values']['dynamic_pressures']*units_dict['pressure']
- densities_values = snap['values']['densities']*units_dict['density']
- masses_values = snap['values']['masses']*units_dict['mass']
+ pressures_values = snap["values"]["pressures"] * units_dict["pressure"]
+ magnetic_pressures_values = (
+ snap["values"]["magnetic_pressures"] * units_dict["pressure"]
+ )
+ dynamic_pressures_values = (
+ snap["values"]["dynamic_pressures"] * units_dict["pressure"]
+ )
+ densities_values = snap["values"]["densities"] * units_dict["density"]
+ masses_values = snap["values"]["masses"] * units_dict["mass"]
Np = len(pressures_values)
if region_cut:
- z = snap['values']['CM_frame_coordinates'][:,2]*units_dict['length']
- y = snap['values']['CM_frame_coordinates'][:,1]*units_dict['length']
- x = snap['values']['CM_frame_coordinates'][:,0]*units_dict['length']
- r = np.sqrt(x**2+y**2)
- mask =((r<=rcut) & (np.abs(z)<=zcut))
+ z = snap["values"]["CM_frame_coordinates"][:, 2] * units_dict["length"]
+ y = snap["values"]["CM_frame_coordinates"][:, 1] * units_dict["length"]
+ x = snap["values"]["CM_frame_coordinates"][:, 0] * units_dict["length"]
+ r = np.sqrt(x ** 2 + y ** 2)
+ mask = (r <= rcut) & (np.abs(z) <= zcut)
pressures_values = pressures_values[mask]
magnetic_pressures_values = magnetic_pressures_values[mask]
dynamic_pressures_values = dynamic_pressures_values[mask]
@@ -405,75 +518,89 @@ def plot_pressure_vs_time(all_snapshots,outputfileaddr,isvec=False,region_cut=Tr
# calculate volume weighted averages:
volume = np.sum(masses_values / densities_values)
- pressures_vv = np.sum(pressures_values * masses_values / densities_values)/volume
- magnetic_pressures_vv = np.sum(magnetic_pressures_values * masses_values / densities_values)/volume
- dynamic_pressures_vv = np.sum(dynamic_pressures_values * masses_values / densities_values)/volume
-
- pressures_vv = pressures_vv.to(units_dict['pressure']).value
- magnetic_pressures_vv = magnetic_pressures_vv.to(units_dict['pressure']).value
- dynamic_pressures_vv = dynamic_pressures_vv.to(units_dict['pressure']).value
+ pressures_vv = (
+ np.sum(pressures_values * masses_values / densities_values) / volume
+ )
+ magnetic_pressures_vv = (
+ np.sum(magnetic_pressures_values * masses_values / densities_values)
+ / volume
+ )
+ dynamic_pressures_vv = (
+ np.sum(dynamic_pressures_values * masses_values / densities_values) / volume
+ )
+
+ pressures_vv = pressures_vv.to(units_dict["pressure"]).value
+ magnetic_pressures_vv = magnetic_pressures_vv.to(units_dict["pressure"]).value
+ dynamic_pressures_vv = dynamic_pressures_vv.to(units_dict["pressure"]).value
pressures_array.append(pressures_vv)
magnetic_pressures_array.append(magnetic_pressures_vv)
dynamic_pressures_array.append(dynamic_pressures_vv)
- time = snap['parameters']['time']
- times_array.append(time)
+ time = snap["parameters"]["time"]
+ times_array.append(time)
npart_array.append(len(pressures_values))
# Retrieve some information about the simulation run
- artDiffusion = snap['metadata']['artDiffusion']
- dedHyp = snap['metadata']['dedHyp']
- dedHypDivv = snap['metadata']['dedHypDivv']
- dedPar = snap['metadata']['dedPar']
- eta = snap['metadata']['eta']
- git = snap['metadata']['git']
- gitBranch = snap['metadata']['gitBranch']
- hydroScheme = snap['metadata']['hydroScheme']
- kernel = snap['metadata']['kernel']
- neighbours = snap['metadata']['neighbours']
-
-
+ artDiffusion = snap["metadata"]["artDiffusion"]
+ dedHyp = snap["metadata"]["dedHyp"]
+ dedHypDivv = snap["metadata"]["dedHypDivv"]
+ dedPar = snap["metadata"]["dedPar"]
+ eta = snap["metadata"]["eta"]
+ git = snap["metadata"]["git"]
+ gitBranch = snap["metadata"]["gitBranch"]
+ hydroScheme = snap["metadata"]["hydroScheme"]
+ kernel = snap["metadata"]["kernel"]
+ neighbours = snap["metadata"]["neighbours"]
+
pressures_array = np.array(pressures_array)
magnetic_pressures_array = np.array(magnetic_pressures_array)
dynamic_pressures_array = np.array(dynamic_pressures_array)
times_array = np.array(times_array)
npart_array = np.array(npart_array)
+ fig, ax = plt.subplots(1, 3, figsize=(5.5 * 3, 5))
+ ax[0].plot(
+ times_array, pressures_array, color="red", marker=".", label=r"$P_{thermal}$"
+ )
+ ax[0].plot(
+ times_array,
+ magnetic_pressures_array,
+ color="blue",
+ marker=".",
+ label=r"$P_{mag}$",
+ )
+ ax[0].plot(
+ times_array,
+ dynamic_pressures_array,
+ color="green",
+ marker=".",
+ label=r"$\rho v^2$",
+ )
- fig, ax = plt.subplots(1, 3, figsize=(5.5*3, 5))
- ax[0].plot(times_array, pressures_array, color='red', marker='.',label=r'$P_{thermal}$')
- ax[0].plot(times_array, magnetic_pressures_array, color='blue', marker='.',label=r'$P_{mag}$')
- ax[0].plot(times_array, dynamic_pressures_array, color='green', marker='.',label=r'$\rho v^2$')
-
ax[0].legend()
ax[0].grid()
- ax[0].set_yscale('log')
- ax[0].set_ylabel('$P$, [$ 10^{10}$ $M_{sol}$ $kpc^{-1}$ $Gyr^{-2}$]',fontsize=20)
- ax[0].set_xlabel('time $[Gyr]$',fontsize=20)
- ax[0].set_ylim([1e-4,1e2])
- ax[0].set_xlim([0,4.0])
- #ax[0].set_ylim([1e-8,1e-4])
-
- ax[1].plot(times_array,npart_array,color='black',marker='.')
+ ax[0].set_yscale("log")
+ ax[0].set_ylabel("$P$, [$ 10^{10}$ $M_{sol}$ $kpc^{-1}$ $Gyr^{-2}$]", fontsize=20)
+ ax[0].set_xlabel("time $[Gyr]$", fontsize=20)
+ ax[0].set_ylim([1e-4, 1e2])
+ ax[0].set_xlim([0, 4.0])
+ # ax[0].set_ylim([1e-8,1e-4])
+
+ ax[1].plot(times_array, npart_array, color="black", marker=".")
ax[1].grid()
- ax[1].set_yscale('log')
- ax[1].set_ylabel(r'$N_p^{cut}$',fontsize=20) #, Lcut = {Lcut} A.U.',fontsize=20)
- ax[1].set_xlabel('time $[Gyr]$',fontsize=20)
- ax[1].set_ylim([1e1,1e5])
- ax[1].set_xlim([0,4.0])
+ ax[1].set_yscale("log")
+ ax[1].set_ylabel(r"$N_p^{cut}$", fontsize=20) # , Lcut = {Lcut} A.U.',fontsize=20)
+ ax[1].set_xlabel("time $[Gyr]$", fontsize=20)
+ ax[1].set_ylim([1e1, 1e5])
+ ax[1].set_xlim([0, 4.0])
# add panel with infromation about the run
text_common_args = dict(
fontsize=10, ha="center", va="center", transform=ax[2].transAxes
)
- ax[2].text(
- 0.5,
- 0.8,
- "Cooling halo with spin, $Np=%.0f$ " % Np,
- **text_common_args,
- )
+ ax[2].text(0.5, 0.8, "Cooling halo with spin, $Np=%.0f$ " % Np, **text_common_args)
ax[2].text(0.5, 0.7, "swift %s" % git.decode("utf-8"), **text_common_args)
ax[2].text(0.5, 0.6, "Branch %s" % gitBranch.decode("utf-8"), **text_common_args)
ax[2].text(0.5, 0.5, hydroScheme.decode("utf-8"), **text_common_args)
@@ -497,111 +624,114 @@ def plot_pressure_vs_time(all_snapshots,outputfileaddr,isvec=False,region_cut=Tr
)
ax[2].axis("off")
-
-
- #fig.suptitle(quantity,fontsize=20)
+ # fig.suptitle(quantity,fontsize=20)
fig.tight_layout()
plt.savefig(outputfileaddr)
plt.close(fig)
- return
+ return
-def plot_B_vs_time(all_snapshots,outputfileaddr,isvec=False,region_cut=True,label='',rcut=15):
+def plot_B_vs_time(
+ all_snapshots, outputfileaddr, isvec=False, region_cut=True, label="", rcut=15
+):
magnetic_flux_densities_array = []
times_array = []
npart_array = []
total_mass_array = []
- rcut*=units_dict['length']
+ rcut *= units_dict["length"]
for snap in all_snapshots:
- magnetic_flux_densities_values = snap['values']['magnetic_flux_densities']*units_dict['magneticfield']
- magnetic_flux_densities_values_sq = dot_vec(magnetic_flux_densities_values,magnetic_flux_densities_values)
- densities_values = snap['values']['densities']*units_dict['density']
- masses_values = snap['values']['masses']*units_dict['mass']
+ magnetic_flux_densities_values = (
+ snap["values"]["magnetic_flux_densities"] * units_dict["magneticfield"]
+ )
+ magnetic_flux_densities_values_sq = dot_vec(
+ magnetic_flux_densities_values, magnetic_flux_densities_values
+ )
+ densities_values = snap["values"]["densities"] * units_dict["density"]
+ masses_values = snap["values"]["masses"] * units_dict["mass"]
Np = len(masses_values)
if region_cut:
- z = snap['values']['CM_frame_coordinates'][:,2]*units_dict['length']
- y = snap['values']['CM_frame_coordinates'][:,1]*units_dict['length']
- x = snap['values']['CM_frame_coordinates'][:,0]*units_dict['length']
- r = np.sqrt(x**2+y**2+z**2)
- mask =(r<=rcut)
+ z = snap["values"]["CM_frame_coordinates"][:, 2] * units_dict["length"]
+ y = snap["values"]["CM_frame_coordinates"][:, 1] * units_dict["length"]
+ x = snap["values"]["CM_frame_coordinates"][:, 0] * units_dict["length"]
+ r = np.sqrt(x ** 2 + y ** 2 + z ** 2)
+ mask = r <= rcut
magnetic_flux_densities_values_sq = magnetic_flux_densities_values_sq[mask]
densities_values = densities_values[mask]
masses_values = masses_values[mask]
# calculate volume weighted averages:
volume = np.sum(masses_values / densities_values)
- magnetic_flux_densities_sq_vv = np.sum(magnetic_flux_densities_values_sq * masses_values / densities_values)/volume
- magnetic_flux_densities_vv = np.sqrt(magnetic_flux_densities_sq_vv)
+ magnetic_flux_densities_sq_vv = (
+ np.sum(magnetic_flux_densities_values_sq * masses_values / densities_values)
+ / volume
+ )
+ magnetic_flux_densities_vv = np.sqrt(magnetic_flux_densities_sq_vv)
- magnetic_flux_densities_vv = magnetic_flux_densities_vv.to(units_dict['magneticfield']).value
+ magnetic_flux_densities_vv = magnetic_flux_densities_vv.to(
+ units_dict["magneticfield"]
+ ).value
magnetic_flux_densities_array.append(magnetic_flux_densities_vv)
- time = snap['parameters']['time']
- times_array.append(time)
+ time = snap["parameters"]["time"]
+ times_array.append(time)
npart_array.append(len(masses_values))
- total_mass_array.append(np.sum(masses_values.to(units_dict['mass']).value))
+ total_mass_array.append(np.sum(masses_values.to(units_dict["mass"]).value))
# Retrieve some information about the simulation run
- artDiffusion = snap['metadata']['artDiffusion']
- dedHyp = snap['metadata']['dedHyp']
- dedHypDivv = snap['metadata']['dedHypDivv']
- dedPar = snap['metadata']['dedPar']
- eta = snap['metadata']['eta']
- git = snap['metadata']['git']
- gitBranch = snap['metadata']['gitBranch']
- hydroScheme = snap['metadata']['hydroScheme']
- kernel = snap['metadata']['kernel']
- neighbours = snap['metadata']['neighbours']
-
-
+ artDiffusion = snap["metadata"]["artDiffusion"]
+ dedHyp = snap["metadata"]["dedHyp"]
+ dedHypDivv = snap["metadata"]["dedHypDivv"]
+ dedPar = snap["metadata"]["dedPar"]
+ eta = snap["metadata"]["eta"]
+ git = snap["metadata"]["git"]
+ gitBranch = snap["metadata"]["gitBranch"]
+ hydroScheme = snap["metadata"]["hydroScheme"]
+ kernel = snap["metadata"]["kernel"]
+ neighbours = snap["metadata"]["neighbours"]
+
magnetic_flux_densities_array = np.array(magnetic_flux_densities_array)
times_array = np.array(times_array)
npart_array = np.array(npart_array)
total_mass_array = np.array(total_mass_array)
- fig, ax = plt.subplots(1, 4, figsize=(5.5*4, 5))
- ax[0].plot(times_array, magnetic_flux_densities_array, color='black', marker='.')
-
- #ax[0].legend()
+ fig, ax = plt.subplots(1, 4, figsize=(5.5 * 4, 5))
+ ax[0].plot(times_array, magnetic_flux_densities_array, color="black", marker=".")
+
+ # ax[0].legend()
ax[0].grid()
- #ax[0].set_yscale('log')
- ax[0].set_ylabel('$B_{rms}$, [$[\mu G]$]',fontsize=20)
- ax[0].set_xlabel('time $[Gyr]$',fontsize=20)
- ax[0].set_ylim([0,3.5])
- ax[0].set_xlim([0,3.0])
- #ax[0].set_ylim([1e-8,1e-4])
-
- ax[1].plot(times_array,npart_array,color='black',marker='.')
+ # ax[0].set_yscale('log')
+ ax[0].set_ylabel("$B_{rms}$, [$[\mu G]$]", fontsize=20)
+ ax[0].set_xlabel("time $[Gyr]$", fontsize=20)
+ ax[0].set_ylim([0, 3.5])
+ ax[0].set_xlim([0, 3.0])
+ # ax[0].set_ylim([1e-8,1e-4])
+
+ ax[1].plot(times_array, npart_array, color="black", marker=".")
ax[1].grid()
- ax[1].set_yscale('log')
- ax[1].set_ylabel(r'$N_p^{cut}$',fontsize=20) #, Lcut = {Lcut} A.U.',fontsize=20)
- ax[1].set_xlabel('time $[Gyr]$',fontsize=20)
- ax[1].set_ylim([1e1,1e5])
- ax[1].set_xlim([0,3.0])
+ ax[1].set_yscale("log")
+ ax[1].set_ylabel(r"$N_p^{cut}$", fontsize=20) # , Lcut = {Lcut} A.U.',fontsize=20)
+ ax[1].set_xlabel("time $[Gyr]$", fontsize=20)
+ ax[1].set_ylim([1e1, 1e5])
+ ax[1].set_xlim([0, 3.0])
-
- ax[2].plot(times_array,total_mass_array,color='black',marker='.')
+ ax[2].plot(times_array, total_mass_array, color="black", marker=".")
ax[2].grid()
- ax[2].set_yscale('log')
- ax[2].set_ylabel(r'$M_{tot} [M_{sol}]$',fontsize=20) #, Lcut = {Lcut} A.U.',fontsize=20)
- ax[2].set_xlabel('time $[Gyr]$',fontsize=20)
- ax[2].set_ylim([1e8,1e11])
- ax[2].set_xlim([0,3.0])
-
+ ax[2].set_yscale("log")
+ ax[2].set_ylabel(
+ r"$M_{tot} [M_{sol}]$", fontsize=20
+ ) # , Lcut = {Lcut} A.U.',fontsize=20)
+ ax[2].set_xlabel("time $[Gyr]$", fontsize=20)
+ ax[2].set_ylim([1e8, 1e11])
+ ax[2].set_xlim([0, 3.0])
# add panel with infromation about the run
text_common_args = dict(
fontsize=10, ha="center", va="center", transform=ax[3].transAxes
)
- ax[3].text(
- 0.5,
- 0.8,
- "Cooling halo with spin, $Np=%.0f$ " % Np,
- **text_common_args,
- )
+ ax[3].text(0.5, 0.8, "Cooling halo with spin, $Np=%.0f$ " % Np, **text_common_args)
ax[3].text(0.5, 0.7, "swift %s" % git.decode("utf-8"), **text_common_args)
ax[3].text(0.5, 0.6, "Branch %s" % gitBranch.decode("utf-8"), **text_common_args)
ax[3].text(0.5, 0.5, hydroScheme.decode("utf-8"), **text_common_args)
@@ -625,23 +755,35 @@ def plot_B_vs_time(all_snapshots,outputfileaddr,isvec=False,region_cut=True,labe
)
ax[3].axis("off")
-
-
- #fig.suptitle(quantity,fontsize=20)
+ # fig.suptitle(quantity,fontsize=20)
fig.tight_layout()
plt.savefig(outputfileaddr)
plt.close(fig)
- return
-
+ return
-# 34013
+# 34013
-folder = './gb099172f_1e6p_hyp=0.1_corr_mass_prof'
+folder = "./gb099172f_1e6p_hyp=0.1_corr_mass_prof"
snapshots_history, snapshot_names = updoad_shapshots_data_from_folder(folder)
-#plot_quantity_vs_time(snapshots_history,'magnetic_flux_densities_norm',outputfileaddr=folder+'/momentums_noRegCut.png',isvec=True,region_cut=False, label = ' All')
-
-plot_pressure_vs_time(snapshots_history,outputfileaddr=folder+'/P_cut.png',isvec=True,region_cut=True,rcut=15,zcut=15.0, label = f' $Lbox=${2*Lcut} kPc')
-plot_B_vs_time(snapshots_history,outputfileaddr=folder+'/B_cut.png',isvec=True,region_cut=True,rcut=15, label = f' $Lbox=${2*Lcut} kPc')
+# plot_quantity_vs_time(snapshots_history,'magnetic_flux_densities_norm',outputfileaddr=folder+'/momentums_noRegCut.png',isvec=True,region_cut=False, label = ' All')
+
+plot_pressure_vs_time(
+ snapshots_history,
+ outputfileaddr=folder + "/P_cut.png",
+ isvec=True,
+ region_cut=True,
+ rcut=15,
+ zcut=15.0,
+ label=f" $Lbox=${2*Lcut} kPc",
+)
+plot_B_vs_time(
+ snapshots_history,
+ outputfileaddr=folder + "/B_cut.png",
+ isvec=True,
+ region_cut=True,
+ rcut=15,
+ label=f" $Lbox=${2*Lcut} kPc",
+)
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/plotCorrelation.py b/examples/MHDTests/CoolingHaloWithSpin_MHD/plotCorrelation.py
index c4199f13afd41dd38ed6910cdd632f93f1914afb..d991ef8a1805fb3be113006ca6c9fc66edeb3e67 100644
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/plotCorrelation.py
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/plotCorrelation.py
@@ -15,9 +15,9 @@ import unyt
filename = sys.argv[1]
-#slice_height = sys.argv[3]
+# slice_height = sys.argv[3]
-#projection = sys.argv[4]
+# projection = sys.argv[4]
prefered_color = "magma"
@@ -46,7 +46,7 @@ y = data.gas.coordinates[:, 1].to_physical()
z = data.gas.coordinates[:, 2].to_physical()
rho = data.gas.densities
-nH = rho.to(unyt.g/unyt.cm**3)/(1.67e-24*unyt.g)
+nH = rho.to(unyt.g / unyt.cm ** 3) / (1.67e-24 * unyt.g)
h = data.gas.smoothing_lengths
v = data.gas.velocities
P = data.gas.pressures
@@ -58,7 +58,7 @@ R1 = data.gas.r1
R2 = data.gas.r2
R3 = data.gas.r3
-normB = np.sqrt(B[:,0]**2+B[:,1]**2+B[:,2]**2)
+normB = np.sqrt(B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2)
# Retrieve some information about the simulation run
artDiffusion = data.metadata.hydro_scheme["Artificial Diffusion Constant"]
@@ -75,21 +75,21 @@ neighbours = data.metadata.hydro_scheme["Kernel target N_ngb"]
########################################## Make histogram
if to_plot == "hist":
- #min_metric = 0.1
+ # min_metric = 0.1
quantity = (
- nH #T.to(unyt.K)#rho.to(unyt.g/unyt.cm**3)
+ nH # T.to(unyt.K)#rho.to(unyt.g/unyt.cm**3)
).value # (rho.to_physical()/rho_mean).value#(abs_vec(v).to_physical()/vrms).value
- qname = r"$n_H$ $[cm^{-3}]$"#r"$\rho$ $[g/cm^3]$" #'rho/<rho>' #'$|v|/v_{rms}$'
+ qname = r"$n_H$ $[cm^{-3}]$" # r"$\rho$ $[g/cm^3]$" #'rho/<rho>' #'$|v|/v_{rms}$'
Nbins = 100
# bins = np.linspace(np.min(quantity),np.max(quantity),Nbins)
bins = np.logspace(
int(np.log10(np.min(quantity))) - 1, int(np.log10(np.max(quantity))) + 1, Nbins
)
fig, ax = plt.subplots()
-
+
nx = 2
ny = 1
- fig, ax = plt.subplots(ny, nx, sharey=True, figsize=(5*nx, 5*ny))
+ fig, ax = plt.subplots(ny, nx, sharey=True, figsize=(5 * nx, 5 * ny))
ax[0].hist(quantity, bins=bins, color="blue", label=qname, alpha=0.5, density=False)
ax[0].set_xscale("log")
@@ -97,11 +97,10 @@ if to_plot == "hist":
ax[0].legend()
ax[0].set_ylabel("$N_{particles}$/bin")
ax[0].set_xlabel(qname)
- ax[0].set_ylim([1,1e5])
- #ax.set_title(f"z={z_to_display:}")
+ ax[0].set_ylim([1, 1e5])
+ # ax.set_title(f"z={z_to_display:}")
ax[0].grid()
-
# add panel with infromation about the run
Np = len(quantity)
text_common_args = dict(
@@ -156,36 +155,35 @@ def prepare_histogram_density_plot(data_x, data_y, Nbins):
)
hist, xedges, yedges = np.histogram2d(data_x, data_y, bins=[bins_x, bins_y])
# building uncorrelated histogram
- #hist_x, _ = np.histogram(data_x, bins=bins_x)
- #hist_y, _ = np.histogram(data_y, bins=bins_y)
- #hist_uncorr = np.outer(hist_x, hist_y).astype("float64")
+ # hist_x, _ = np.histogram(data_x, bins=bins_x)
+ # hist_y, _ = np.histogram(data_y, bins=bins_y)
+ # hist_uncorr = np.outer(hist_x, hist_y).astype("float64")
# norm
min_level = 1 / np.sum(hist)
hist /= np.sum(hist)
- #hist_uncorr /= np.sum(hist_uncorr)
+ # hist_uncorr /= np.sum(hist_uncorr)
# print(np.sum(hist),np.sum(hist_uncorr))
# regularize
- #hist[hist < min_level] = min_level
- #hist_uncorr[hist_uncorr < min_level] = min_level
+ # hist[hist < min_level] = min_level
+ # hist_uncorr[hist_uncorr < min_level] = min_level
# build correlation_funciton
- #correlation_funcion = hist - hist_uncorr
+ # correlation_funcion = hist - hist_uncorr
X, Y = np.meshgrid(xedges, yedges)
- return X,Y, hist # hist_uncorr, correlation_funcion
+ return X, Y, hist # hist_uncorr, correlation_funcion
+
def plot_histogram_density_plot(quantity_x, quantity_y, qname_x, qname_y, Nbins=100):
quantity_x = quantity_x.value
quantity_y = quantity_y.value
-
- X,Y,hist = prepare_histogram_density_plot(
- quantity_x, quantity_y, Nbins
- )
+
+ X, Y, hist = prepare_histogram_density_plot(quantity_x, quantity_y, Nbins)
nx = 2
ny = 1
- fig, ax = plt.subplots(ny, nx, sharey=True, figsize=(5*nx, 5*ny))
+ fig, ax = plt.subplots(ny, nx, sharey=True, figsize=(5 * nx, 5 * ny))
fig.suptitle(f"({qname_x},{qname_y}) space")
@@ -193,17 +191,17 @@ def plot_histogram_density_plot(quantity_x, quantity_y, qname_x, qname_y, Nbins=
ax[0].set_ylabel(qname_y)
ax[0].set_xlabel(qname_x)
ax[0].set_box_aspect(1)
- #plt.imshow(hist.T, origin='lower', aspect='auto', cmap='plasma')
+ # plt.imshow(hist.T, origin='lower', aspect='auto', cmap='plasma')
pax = ax[0].pcolormesh(X, Y, hist.T, cmap=prefered_color, norm=LogNorm())
- ax[0].grid(True, linewidth=0.5,alpha=0.3)
+ ax[0].grid(True, linewidth=0.5, alpha=0.3)
ax[0].set_xscale("log")
ax[0].set_yscale("log")
- x = np.logspace(np.log10(min(quantity_x)),np.log10(max(quantity_x)),100)
- y = min(quantity_y)/10 * (x/min(quantity_x))**(2/3)
+ x = np.logspace(np.log10(min(quantity_x)), np.log10(max(quantity_x)), 100)
+ y = min(quantity_y) / 10 * (x / min(quantity_x)) ** (2 / 3)
- #ax[0].plot(x,y,color='red',alpha=0.5)
+ # ax[0].plot(x,y,color='red',alpha=0.5)
# add panel with infromation about the run
Np = len(quantity_x)
@@ -243,17 +241,19 @@ def plot_histogram_density_plot(quantity_x, quantity_y, qname_x, qname_y, Nbins=
)
ax[1].axis("off")
- #fig.colorbar(pax)
+ # fig.colorbar(pax)
fig.tight_layout()
plt.savefig(sys.argv[2], dpi=200)
+
if to_plot == "correlation_hist":
- q_x = nH.to(1/unyt.cm**3)
- q_y = T.to(unyt.K) #normB.to(1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s)) #T.to(unyt.K)
+ q_x = nH.to(1 / unyt.cm ** 3)
+ q_y = T.to(
+ unyt.K
+ ) # normB.to(1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s)) #T.to(unyt.K)
q_x = q_x.to_physical()
q_y = q_y.to_physical()
plot_histogram_density_plot(
q_x, q_y, r"$\rho / m_H$ $[1/cm^3]$", "T $[K]$", Nbins=200
)
exit()
-
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/plotErrors.py b/examples/MHDTests/CoolingHaloWithSpin_MHD/plotErrors.py
index a9a56483b183ac763f94d23732b58f53589087fe..bbb6d9d639245c854eb0ab3aa0f4b7d22a32e729 100644
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/plotErrors.py
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/plotErrors.py
@@ -40,13 +40,13 @@ v = data.gas.velocities
norm_r = np.sqrt(r[:, 0] ** 2 + r[:, 1] ** 2 + r[:, 2] ** 2)
vr = np.sum(v * r, axis=1) / norm_r
normv = np.sqrt(v[:, 0] ** 2 + v[:, 1] ** 2 + v[:, 2] ** 2)
-B = data.gas.magnetic_flux_densities
-#B = v.value * 0.0 * 1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s)
-#J = data.gas.magnetic_flux_curl
+B = data.gas.magnetic_flux_densities
+# B = v.value * 0.0 * 1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s)
+# J = data.gas.magnetic_flux_curl
normB = np.sqrt(B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2)
-#normJ = np.sqrt(J[:, 0] ** 2 + J[:, 1] ** 2 + J[:, 2] ** 2)
+# normJ = np.sqrt(J[:, 0] ** 2 + J[:, 1] ** 2 + J[:, 2] ** 2)
divB = data.gas.magnetic_divergences
-#divB = norm_r.value * 0.0 * 1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s* unyt.cm)
+# divB = norm_r.value * 0.0 * 1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s* unyt.cm)
h = data.gas.smoothing_lengths
rho = data.gas.densities
@@ -57,22 +57,21 @@ R2 = data.gas.r2
OW = data.gas.owtriggers
mean_OW = np.mean(OW)
min_OW = np.min(OW)
-print('Mean OW is ',np.round(mean_OW,1))
-print('Min OW is ', np.round(min_OW,1))
-OW = OW/mean_OW
+print("Mean OW is ", np.round(mean_OW, 1))
+print("Min OW is ", np.round(min_OW, 1))
+OW = OW / mean_OW
# Normalize errors
-R0[R0<1e-2] = 1e-2
-R1[R1<1e-2] = 1e-2
-R2[R2<1e-2] = 1e-2
-OW[OW<1e-2] = 1e-2
-
+R0[R0 < 1e-2] = 1e-2
+R1[R1 < 1e-2] = 1e-2
+R2[R2 < 1e-2] = 1e-2
+OW[OW < 1e-2] = 1e-2
Np = len(h)
-print('Number of particles %E' % len(data.gas.masses))
-print('Total mass %E Msol' % np.sum(data.gas.masses.to(unyt.Msun).value))
-print('Gas particle mass %E Msol' % np.mean(data.gas.masses.to(unyt.Msun).value))
+print("Number of particles %E" % len(data.gas.masses))
+print("Total mass %E Msol" % np.sum(data.gas.masses.to(unyt.Msun).value))
+print("Gas particle mass %E Msol" % np.mean(data.gas.masses.to(unyt.Msun).value))
data.gas.mass_weighted_R0 = data.gas.masses * R0
data.gas.mass_weighted_R1 = data.gas.masses * R1
@@ -88,7 +87,7 @@ rotation_matrix = [[1,0,0],
rotation_matrix = [[0, 0, 1], [0, 1, 0], [-1, 0, 0]]
# set plot area
-Lslice_kPc = 20.0 #2*21.5
+Lslice_kPc = 20.0 # 2*21.5
Lslice = Lslice_kPc * 3.08e18 * 1e3 * unyt.cm
visualise_region_xy = [
@@ -139,7 +138,7 @@ mass_weighted_R1_map = slice_gas(**common_arguments, project="mass_weighted_R1")
mass_weighted_R2_map = slice_gas(**common_arguments, project="mass_weighted_R2")
mass_weighted_OW_map = slice_gas(**common_arguments, project="mass_weighted_OW")
-#mass_weighted_J_map = slice_gas(**common_arguments, project="mass_weighted_J")
+# mass_weighted_J_map = slice_gas(**common_arguments, project="mass_weighted_J")
R0_map_xy = mass_weighted_R0_map_xy / mass_map_xy
R1_map_xy = mass_weighted_R1_map_xy / mass_map_xy
R2_map_xy = mass_weighted_R2_map_xy / mass_map_xy
@@ -152,60 +151,60 @@ OW_map = mass_weighted_OW_map / mass_map
nx = 2
ny = 5
-fig, ax = plt.subplots(ny, nx, sharey=True, figsize=(5*nx, 5*ny))
+fig, ax = plt.subplots(ny, nx, sharey=True, figsize=(5 * nx, 5 * ny))
-#fig.suptitle("Plotting at %.3f free fall times" % toplot)
+# fig.suptitle("Plotting at %.3f free fall times" % toplot)
a00 = ax[0, 0].contourf(
np.log10(R0_map.value),
- cmap="plasma",#"plasma",#"gist_stern",#"RdBu",
+ cmap="plasma", # "plasma",#"gist_stern",#"RdBu",
extend="both",
levels=np.linspace(-2, 0, 100),
)
a01 = ax[0, 1].contourf(
np.log10(R0_map_xy.value).T,
- cmap="plasma",#"plasma",#"gist_stern", #"RdBu",
+ cmap="plasma", # "plasma",#"gist_stern", #"RdBu",
extend="both",
levels=np.linspace(-2, 0, 100),
)
a10 = ax[1, 0].contourf(
np.log10(R1_map.value),
- cmap="plasma",#"viridis", #"nipy_spectral_r",
+ cmap="plasma", # "viridis", #"nipy_spectral_r",
extend="both",
levels=np.linspace(-2, 0, 100),
)
a11 = ax[1, 1].contourf(
np.log10(R1_map_xy.value).T,
- cmap="plasma",#"viridis", #"nipy_spectral_r",
+ cmap="plasma", # "viridis", #"nipy_spectral_r",
extend="both",
levels=np.linspace(-2, 0, 100),
)
a20 = ax[2, 0].contourf(
np.log10(R2_map.value),
- cmap="plasma",#"viridis", #"nipy_spectral_r",
+ cmap="plasma", # "viridis", #"nipy_spectral_r",
extend="both",
levels=np.linspace(-2, 0, 100),
)
a21 = ax[2, 1].contourf(
np.log10(R2_map_xy.value).T,
- cmap="plasma",#"viridis", #"nipy_spectral_r",
+ cmap="plasma", # "viridis", #"nipy_spectral_r",
extend="both",
levels=np.linspace(-2, 0, 100),
)
a30 = ax[3, 0].contourf(
np.log10(OW_map.value),
- cmap="bwr",#"viridis", #"nipy_spectral_r",
+ cmap="bwr", # "viridis", #"nipy_spectral_r",
extend="both",
levels=np.linspace(-1, 1, 100),
)
a31 = ax[3, 1].contourf(
np.log10(OW_map_xy.value).T,
- cmap="bwr",#"viridis", #"nipy_spectral_r",
+ cmap="bwr", # "viridis", #"nipy_spectral_r",
extend="both",
levels=np.linspace(-1, 1, 100),
)
@@ -222,45 +221,37 @@ for ii in range(ny):
for jj in range(nx):
ax[ii, jj].set_xlabel(r"$x$ [kPc]")
ax[ii, jj].set_xticks(locs, labels)
- #ax[ii, jj].set_aspect("equal")
- ax[ii,jj].set_aspect('equal', adjustable='box')
- ax[ii,jj].set_xlim(0, 511)
- ax[ii,jj].set_ylim(0, 511)
+ # ax[ii, jj].set_aspect("equal")
+ ax[ii, jj].set_aspect("equal", adjustable="box")
+ ax[ii, jj].set_xlim(0, 511)
+ ax[ii, jj].set_ylim(0, 511)
-ticks_error = [-2,-1,0]
-cbar1 = plt.colorbar(
- a00, ax=ax[0, 0], fraction=0.046, pad=0.04, ticks=ticks_error
-)
+ticks_error = [-2, -1, 0]
+cbar1 = plt.colorbar(a00, ax=ax[0, 0], fraction=0.046, pad=0.04, ticks=ticks_error)
cbar1.set_label(r"$\mathrm{log}_{10} R_0 \; $")
-cbar2 = plt.colorbar(
- a01, ax=ax[0, 1], fraction=0.046, pad=0.04, ticks=ticks_error
-)
+cbar2 = plt.colorbar(a01, ax=ax[0, 1], fraction=0.046, pad=0.04, ticks=ticks_error)
ax[0, 1].set_ylabel(r"$y$ [kPc]")
cbar2.set_label(r"$\mathrm{log}_{10} R_0 \; $")
ax[1, 1].set_ylabel(r"$y$ [kPc]")
-cbar3 = plt.colorbar(
- a11, ax=ax[1, 0], fraction=0.046, pad=0.04, ticks=ticks_error
-)
+cbar3 = plt.colorbar(a11, ax=ax[1, 0], fraction=0.046, pad=0.04, ticks=ticks_error)
cbar3.set_label(r"$\mathrm{log}_{10} R_1 \;$")
-cbar4 = plt.colorbar(
- a11, ax=ax[1, 1], fraction=0.046, pad=0.04, ticks=ticks_error
-)
+cbar4 = plt.colorbar(a11, ax=ax[1, 1], fraction=0.046, pad=0.04, ticks=ticks_error)
cbar4.set_label(r"$\mathrm{log}_{10} R_1 \;$")
ax[2, 1].set_ylabel(r"$y$ [kPc]")
-cbar5 = plt.colorbar(a20, ax=ax[2, 0], fraction=0.046, pad=0.04, ticks=ticks_error)
+cbar5 = plt.colorbar(a20, ax=ax[2, 0], fraction=0.046, pad=0.04, ticks=ticks_error)
cbar5.set_label(r"$\mathrm{log}_{10} R_2 \;$")
cbar6 = plt.colorbar(a21, ax=ax[2, 1], fraction=0.046, pad=0.04, ticks=ticks_error)
cbar6.set_label(r"$\mathrm{log}_{10} R_2 \;$")
-ticksOW=[-1,0,1]
-cbar7 = plt.colorbar(a30, ax=ax[3, 0], fraction=0.046, pad=0.04,ticks=ticksOW)
+ticksOW = [-1, 0, 1]
+cbar7 = plt.colorbar(a30, ax=ax[3, 0], fraction=0.046, pad=0.04, ticks=ticksOW)
cbar7.set_label(r"$\mathrm{log}_{10} OW/\langle OW \rangle \;$")
-cbar8 = plt.colorbar(a31, ax=ax[3, 1], fraction=0.046, pad=0.04,ticks=ticksOW)
+cbar8 = plt.colorbar(a31, ax=ax[3, 1], fraction=0.046, pad=0.04, ticks=ticksOW)
cbar8.set_label(r"$\mathrm{log}_{10} OW/\langle OW \rangle \;$")
Ninfo = 4
@@ -298,11 +289,14 @@ ax[Ninfo, 1].text(
0.5, 0.1, "Physical resistivity $\eta$: $%.2f$ " % (eta), **text_common_args
)
ax[Ninfo, 1].text(
- 0.5, 0.0, r"Number of particles $N_p$: $%.0f$, $\langle OW \rangle$ = $%.0f$, $ OW_{min}$ = $%.0f$ " % (Np,mean_OW,min_OW), **text_common_args
+ 0.5,
+ 0.0,
+ r"Number of particles $N_p$: $%.0f$, $\langle OW \rangle$ = $%.0f$, $ OW_{min}$ = $%.0f$ "
+ % (Np, mean_OW, min_OW),
+ **text_common_args,
)
ax[Ninfo, 1].axis("off")
fig.tight_layout()
-plt.savefig(sys.argv[2],dpi=220)
-
+plt.savefig(sys.argv[2], dpi=220)
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/plotMovie.py b/examples/MHDTests/CoolingHaloWithSpin_MHD/plotMovie.py
index 7126b177cf3242e6455da9cfddf2060e33ac634f..a364c6f574230fda3a8e9f10fdd925c0b4e9dc60 100644
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/plotMovie.py
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/plotMovie.py
@@ -1,30 +1,36 @@
import subprocess
from glob import glob
-import os
+import os
+
def import_all_snapshot_addr(folderaddr):
- addr_book = glob(folderaddr+'**/CoolingHalo_*.hdf5',recursive=True)
+ addr_book = glob(folderaddr + "**/CoolingHalo_*.hdf5", recursive=True)
addr_book = sorted(addr_book)
- image_names = ['image_'+addr.split('.hdf5')[0].split('_')[-1]+'.png' for addr in addr_book]
+ image_names = [
+ "image_" + addr.split(".hdf5")[0].split("_")[-1] + ".png" for addr in addr_book
+ ]
return addr_book, image_names
+
def plot_all_snapshots(folderaddr):
snap_addr, output_addr = import_all_snapshot_addr(folderaddr)
for saddr, iaddr in zip(snap_addr, output_addr):
- command = 'python3 plotCorrelation.py '+ saddr + ' ' + iaddr
+ command = "python3 plotCorrelation.py " + saddr + " " + iaddr
try:
- subprocess.call(command,shell=True)
+ subprocess.call(command, shell=True)
except subprocess.CalledProcessError as e:
- print(f'Encountered error number {e.returncode}, command output: {e.output}')
+ print(
+ f"Encountered error number {e.returncode}, command output: {e.output}"
+ )
-def make_a_movie(framerate = 1):
- command = f'ffmpeg -framerate {framerate} -i \'image_%04d.png\' -avoid_negative_ts make_zero -r {framerate} -pix_fmt yuv420p video.mp4'
+
+def make_a_movie(framerate=1):
+ command = f"ffmpeg -framerate {framerate} -i 'image_%04d.png' -avoid_negative_ts make_zero -r {framerate} -pix_fmt yuv420p video.mp4"
try:
- subprocess.call(command,shell=True)
+ subprocess.call(command, shell=True)
except subprocess.CalledProcessError as e:
- print(f'Encountered error number {e.returncode}, command output: {e.output}')
-
+ print(f"Encountered error number {e.returncode}, command output: {e.output}")
-plot_all_snapshots('./CH_gb099172f_hyp=0.1_par=1.0_divv=0.0_smallHalo_1e4_correcteps')
-#make_a_movie()
+plot_all_snapshots("./CH_gb099172f_hyp=0.1_par=1.0_divv=0.0_smallHalo_1e4_correcteps")
+# make_a_movie()
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/plotSolutionReference.py b/examples/MHDTests/CoolingHaloWithSpin_MHD/plotSolutionReference.py
index 4938b9eb3842e45ace708027da7b9e0e7fbf36f1..a563d432517d48f16a5a76f48af9729ac6138c01 100644
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/plotSolutionReference.py
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/plotSolutionReference.py
@@ -40,21 +40,21 @@ v = data.gas.velocities
norm_r = np.sqrt(r[:, 0] ** 2 + r[:, 1] ** 2 + r[:, 2] ** 2)
vr = np.sum(v * r, axis=1) / norm_r
normv = np.sqrt(v[:, 0] ** 2 + v[:, 1] ** 2 + v[:, 2] ** 2)
-B = data.gas.magnetic_flux_densities
-#B = v.value * 0.0 * 1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s)
-#J = data.gas.magnetic_flux_curl
+B = data.gas.magnetic_flux_densities
+# B = v.value * 0.0 * 1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s)
+# J = data.gas.magnetic_flux_curl
normB = np.sqrt(B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2)
-#normJ = np.sqrt(J[:, 0] ** 2 + J[:, 1] ** 2 + J[:, 2] ** 2)
+# normJ = np.sqrt(J[:, 0] ** 2 + J[:, 1] ** 2 + J[:, 2] ** 2)
divB = data.gas.magnetic_divergences
-#divB = norm_r.value * 0.0 * 1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s* unyt.cm)
+# divB = norm_r.value * 0.0 * 1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s* unyt.cm)
h = data.gas.smoothing_lengths
rho = data.gas.densities
Np = len(h)
-print('Number of particles %E' % len(data.gas.masses))
-print('Total mass %E Msol' % np.sum(data.gas.masses.to(unyt.Msun).value))
-print('Gas particle mass %E Msol' % np.mean(data.gas.masses.to(unyt.Msun).value))
+print("Number of particles %E" % len(data.gas.masses))
+print("Total mass %E Msol" % np.sum(data.gas.masses.to(unyt.Msun).value))
+print("Gas particle mass %E Msol" % np.mean(data.gas.masses.to(unyt.Msun).value))
-plasmaBeta = data.gas.pressures/(normB**2/(2*mu0))
+plasmaBeta = data.gas.pressures / (normB ** 2 / (2 * mu0))
data.gas.mass_weighted_densities = data.gas.masses * data.gas.densities
data.gas.mass_weighted_error = data.gas.masses * np.log10(
@@ -64,11 +64,11 @@ data.gas.mass_weighted_vr = data.gas.masses * vr
data.gas.mass_weighted_Bz = data.gas.masses * abs(B[:, 2])
data.gas.mass_weighted_Bx = data.gas.masses * abs(B[:, 0])
data.gas.mass_weighted_By = data.gas.masses * abs(B[:, 1])
-#data.gas.mass_weighted_J = (
+# data.gas.mass_weighted_J = (
# data.gas.masses * np.sqrt(J[:, 0] ** 2 + J[:, 1] ** 2 + J[:, 2] ** 2) / mu0
-#)
+# )
-data.gas.mass_weighted_plasmaBeta = data.gas.masses * plasmaBeta
+data.gas.mass_weighted_plasmaBeta = data.gas.masses * plasmaBeta
data.gas.mass_weighted_normv = data.gas.masses * normv
data.gas.mass_weighted_normB = data.gas.masses * normB
@@ -80,7 +80,7 @@ rotation_matrix = [[1,0,0],
rotation_matrix = [[0, 0, 1], [0, 1, 0], [-1, 0, 0]]
# set plot area
-Lslice_kPc = 20.0 #2*21.5
+Lslice_kPc = 20.0 # 2*21.5
Lslice = Lslice_kPc * 3.08e18 * 1e3 * unyt.cm
visualise_region_xy = [
@@ -129,18 +129,26 @@ mass_weighted_density_map = slice_gas(
**common_arguments, project="mass_weighted_densities"
)
mass_weighted_error_map = slice_gas(**common_arguments, project="mass_weighted_error")
-mass_weighted_error_map_xy = slice_gas(**common_arguments_xy, project="mass_weighted_error")
+mass_weighted_error_map_xy = slice_gas(
+ **common_arguments_xy, project="mass_weighted_error"
+)
mass_weighted_vr_map = slice_gas(**common_arguments, project="mass_weighted_vr")
mass_weighted_Bz_map = slice_gas(**common_arguments, project="mass_weighted_Bz")
mass_weighted_Bx_map = slice_gas(**common_arguments, project="mass_weighted_Bx")
mass_weighted_By_map = slice_gas(**common_arguments, project="mass_weighted_By")
-#mass_weighted_J_map = slice_gas(**common_arguments, project="mass_weighted_J")
+# mass_weighted_J_map = slice_gas(**common_arguments, project="mass_weighted_J")
-mass_weighted_plasmaBeta_map = slice_gas(**common_arguments, project="mass_weighted_plasmaBeta")
+mass_weighted_plasmaBeta_map = slice_gas(
+ **common_arguments, project="mass_weighted_plasmaBeta"
+)
mass_weighted_normv_map = slice_gas(**common_arguments, project="mass_weighted_normv")
-mass_weighted_normv_map_xy = slice_gas(**common_arguments_xy, project="mass_weighted_normv")
+mass_weighted_normv_map_xy = slice_gas(
+ **common_arguments_xy, project="mass_weighted_normv"
+)
mass_weighted_normB_map = slice_gas(**common_arguments, project="mass_weighted_normB")
-mass_weighted_normB_map_xy = slice_gas(**common_arguments_xy, project="mass_weighted_normB")
+mass_weighted_normB_map_xy = slice_gas(
+ **common_arguments_xy, project="mass_weighted_normB"
+)
mass_weighted_Bz_map_xy = slice_gas(**common_arguments_xy, project="mass_weighted_Bz")
mass_weighted_Bx_map_xy = slice_gas(**common_arguments_xy, project="mass_weighted_Bx")
mass_weighted_By_map_xy = slice_gas(**common_arguments_xy, project="mass_weighted_By")
@@ -156,85 +164,85 @@ By_map = mass_weighted_By_map / mass_map
Bz_map_xy = mass_weighted_Bz_map_xy / mass_map_xy
Bx_map_xy = mass_weighted_Bx_map_xy / mass_map_xy
By_map_xy = mass_weighted_By_map_xy / mass_map_xy
-#J_map = mass_weighted_J_map / mass_map
+# J_map = mass_weighted_J_map / mass_map
plasmaBeta_map = mass_weighted_plasmaBeta_map / mass_map
normv_map = mass_weighted_normv_map / mass_map
normv_map_xy = mass_weighted_normv_map_xy / mass_map_xy
normB_map = mass_weighted_normB_map / mass_map
normB_map_xy = mass_weighted_normB_map_xy / mass_map_xy
-#density_map_xy.convert_to_units(unyt.g * unyt.cm ** (-3))
-#density_map.convert_to_units(unyt.g * unyt.cm ** (-3))
+# density_map_xy.convert_to_units(unyt.g * unyt.cm ** (-3))
+# density_map.convert_to_units(unyt.g * unyt.cm ** (-3))
-density_map_xy.convert_to_units(1e10*(1e33*unyt.g) * (3.08e21*unyt.cm) ** (-3))
-density_map.convert_to_units(1e10*(1e33*unyt.g) * (3.08e21*unyt.cm) ** (-3))
+density_map_xy.convert_to_units(1e10 * (1e33 * unyt.g) * (3.08e21 * unyt.cm) ** (-3))
+density_map.convert_to_units(1e10 * (1e33 * unyt.g) * (3.08e21 * unyt.cm) ** (-3))
-nH_map_xy = density_map_xy/(1.67e-24 * unyt.g)
-nH_map = density_map/(1.67e-24 * unyt.g)
+nH_map_xy = density_map_xy / (1.67e-24 * unyt.g)
+nH_map = density_map / (1.67e-24 * unyt.g)
vr_map.convert_to_units(unyt.km / unyt.s)
-Bz_map.convert_to_units(1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s))
-Bx_map.convert_to_units(1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s))
-By_map.convert_to_units(1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s))
-Bz_map_xy.convert_to_units(1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s))
-Bx_map_xy.convert_to_units(1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s))
-By_map_xy.convert_to_units(1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s))
-#J_map.convert_to_units(unyt.statA / (unyt.m ** 2))
+Bz_map.convert_to_units(1e-7 * unyt.g / (unyt.statA * unyt.s * unyt.s))
+Bx_map.convert_to_units(1e-7 * unyt.g / (unyt.statA * unyt.s * unyt.s))
+By_map.convert_to_units(1e-7 * unyt.g / (unyt.statA * unyt.s * unyt.s))
+Bz_map_xy.convert_to_units(1e-7 * unyt.g / (unyt.statA * unyt.s * unyt.s))
+Bx_map_xy.convert_to_units(1e-7 * unyt.g / (unyt.statA * unyt.s * unyt.s))
+By_map_xy.convert_to_units(1e-7 * unyt.g / (unyt.statA * unyt.s * unyt.s))
+# J_map.convert_to_units(unyt.statA / (unyt.m ** 2))
normv_map.convert_to_units(unyt.km / unyt.s)
normv_map_xy.convert_to_units(unyt.km / unyt.s)
-normB_map.convert_to_units(1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s))
-normB_map_xy.convert_to_units(1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s))
-#normB_map = np.sqrt(Bx_map**2+By_map**2+Bz_map**2)
-#normB_map_xy = np.sqrt(Bx_map_xy**2+By_map_xy**2+Bz_map_xy**2)
+normB_map.convert_to_units(1e-7 * unyt.g / (unyt.statA * unyt.s * unyt.s))
+normB_map_xy.convert_to_units(1e-7 * unyt.g / (unyt.statA * unyt.s * unyt.s))
+# normB_map = np.sqrt(Bx_map**2+By_map**2+Bz_map**2)
+# normB_map_xy = np.sqrt(Bx_map_xy**2+By_map_xy**2+Bz_map_xy**2)
nx = 2
ny = 4
-fig, ax = plt.subplots(ny, nx, sharey=True, figsize=(5*nx, 5*ny))
+fig, ax = plt.subplots(ny, nx, sharey=True, figsize=(5 * nx, 5 * ny))
-#fig.suptitle("Plotting at %.3f free fall times" % toplot)
+# fig.suptitle("Plotting at %.3f free fall times" % toplot)
a00 = ax[0, 0].contourf(
np.log10(nH_map.value),
- #np.log10(density_map.value),
- cmap="jet",#"plasma",#"gist_stern",#"RdBu",
+ # np.log10(density_map.value),
+ cmap="jet", # "plasma",#"gist_stern",#"RdBu",
extend="both",
levels=np.linspace(-4, 2, 100),
- #levels=np.linspace(-7, -1, 100),
+ # levels=np.linspace(-7, -1, 100),
)
a01 = ax[0, 1].contourf(
np.log10(nH_map_xy.value).T,
- #np.log10(density_map_xy.value).T,
- cmap="jet",#"plasma",#"gist_stern", #"RdBu",
+ # np.log10(density_map_xy.value).T,
+ cmap="jet", # "plasma",#"gist_stern", #"RdBu",
extend="both",
levels=np.linspace(-4, 2, 100),
- #levels=np.linspace(-7, -1, 100),
+ # levels=np.linspace(-7, -1, 100),
)
-#a10 = ax[1, 0].contourf(
+# a10 = ax[1, 0].contourf(
# np.log10(np.maximum(normv_map.value,-1)), cmap="cividis", extend="both", levels=np.linspace(1.0, 2.0, 100)
-#)
-#a11 = ax[1, 1].contourf(
+# )
+# a11 = ax[1, 1].contourf(
# np.log10(np.maximum(normv_map_xy.value,-1)), cmap="cividis", extend="both", levels=np.linspace(1.0, 2.0, 100)
-#)
-#a11 = ax[1, 1].contourf(
+# )
+# a11 = ax[1, 1].contourf(
# error_map.value, cmap="jet", extend="both", levels=np.arange(-6.0, 3.0, 1.0)
-#)
-#a20 = ax[2, 0].contourf(
+# )
+# a20 = ax[2, 0].contourf(
# np.log10(np.maximum(plasmaBeta_map.value, -4)),
# cmap="gist_stern",
# extend="both",
# levels=np.linspace(-3, 4, 100),
-#)
+# )
a10 = ax[1, 0].contourf(
np.maximum(np.log10(normB_map.value), -6),
- cmap="jet",#"viridis", #"nipy_spectral_r",
+ cmap="jet", # "viridis", #"nipy_spectral_r",
extend="both",
levels=np.linspace(-2, 2, 100),
)
a11 = ax[1, 1].contourf(
np.maximum(np.log10(normB_map_xy.value), -6).T,
- cmap="jet",#"viridis", #"nipy_spectral_r",
+ cmap="jet", # "viridis", #"nipy_spectral_r",
extend="both",
levels=np.linspace(-2, 2, 100),
)
@@ -247,7 +255,10 @@ a21 = ax[2, 1].contourf(
)
a30 = ax[3, 0].contourf(
- np.maximum(np.log10(plasmaBeta_map.value),-10), cmap="bwr", extend="both", levels=np.linspace(0.0, 1.0, 100)
+ np.maximum(np.log10(plasmaBeta_map.value), -10),
+ cmap="bwr",
+ extend="both",
+ levels=np.linspace(0.0, 1.0, 100),
)
locs = [512 / 4, 512 / 2, 3 * 512 / 4]
locs = [512 / 4, 512 / 2, 3 * 512 / 4]
@@ -261,51 +272,43 @@ for ii in range(ny):
for jj in range(nx):
ax[ii, jj].set_xlabel(r"$x$ [kPc]")
ax[ii, jj].set_xticks(locs, labels)
- #ax[ii, jj].set_aspect("equal")
- ax[ii,jj].set_aspect('equal', adjustable='box')
- ax[ii,jj].set_xlim(0, 511)
- ax[ii,jj].set_ylim(0, 511)
+ # ax[ii, jj].set_aspect("equal")
+ ax[ii, jj].set_aspect("equal", adjustable="box")
+ ax[ii, jj].set_xlim(0, 511)
+ ax[ii, jj].set_ylim(0, 511)
-ticks_rho = [-4,-3,-2,-1, 0, 1, 2]
-#ticks_rho = [-7,-6,-5,-4,-3,-2,-1]
-cbar1 = plt.colorbar(
- a00, ax=ax[0, 0], fraction=0.046, pad=0.04, ticks=ticks_rho
-)
+ticks_rho = [-4, -3, -2, -1, 0, 1, 2]
+# ticks_rho = [-7,-6,-5,-4,-3,-2,-1]
+cbar1 = plt.colorbar(a00, ax=ax[0, 0], fraction=0.046, pad=0.04, ticks=ticks_rho)
cbar1.set_label(r"$\mathrm{log}_{10} \rho / m_H \; [ {cm}^{-3}]$")
-#cbar1.set_label(r"$\mathrm{log}_{10} \rho \; [10^{10} {M}_{sol} {kpc}^{-3}]$")
-cbar2 = plt.colorbar(
- a01, ax=ax[0, 1], fraction=0.046, pad=0.04, ticks=ticks_rho
-)
+# cbar1.set_label(r"$\mathrm{log}_{10} \rho \; [10^{10} {M}_{sol} {kpc}^{-3}]$")
+cbar2 = plt.colorbar(a01, ax=ax[0, 1], fraction=0.046, pad=0.04, ticks=ticks_rho)
ax[0, 1].set_ylabel(r"$y$ [kPc]")
cbar2.set_label(r"$\mathrm{log}_{10} \rho / m_H \; [ {cm}^{-3}]$")
-#cbar2.set_label(r"$\mathrm{log}_{10} \rho \; [10^{10} {M}_{sol} {kpc}^{-3}]$")
+# cbar2.set_label(r"$\mathrm{log}_{10} \rho \; [10^{10} {M}_{sol} {kpc}^{-3}]$")
-#cbar3 = plt.colorbar(
+# cbar3 = plt.colorbar(
# a10, ax=ax[1, 0], fraction=0.046, pad=0.04, ticks=np.linspace(1, 2, 2)
-#)
-#cbar3.set_label(r"$\mathrm{log}_{10} |v| \; [km s^{-1}]$")
+# )
+# cbar3.set_label(r"$\mathrm{log}_{10} |v| \; [km s^{-1}]$")
#
-#cbar4 = plt.colorbar(
+# cbar4 = plt.colorbar(
# a11, ax=ax[1, 1], fraction=0.046, pad=0.04, ticks=np.linspace(1, 2, 2)
-#)
-#cbar4.set_label(r"$\mathrm{log}_{10} |v| \; [km s^{-1}]$")
+# )
+# cbar4.set_label(r"$\mathrm{log}_{10} |v| \; [km s^{-1}]$")
ax[1, 1].set_ylabel(r"$y$ [kPc]")
-#cbar4 = plt.colorbar(a11, ax=ax[1, 1], fraction=0.046, pad=0.04)
-#cbar4.set_label(r"$\mathrm{log}_{10} \frac{h \: \nabla \cdot B}{|B|}$")
-#ticks_Beta = [-3,-1.83,-0.667,0.500,1.67,2.83,4.0]
-#cbar5 = plt.colorbar(
+# cbar4 = plt.colorbar(a11, ax=ax[1, 1], fraction=0.046, pad=0.04)
+# cbar4.set_label(r"$\mathrm{log}_{10} \frac{h \: \nabla \cdot B}{|B|}$")
+# ticks_Beta = [-3,-1.83,-0.667,0.500,1.67,2.83,4.0]
+# cbar5 = plt.colorbar(
# a20, ax=ax[2, 0], fraction=0.046, pad=0.04, ticks=ticks_Beta
-#)
-#cbar5.set_label(r"$\mathrm{log}_{10} \beta \;$")
-ticks_B = [-2,-1,0,1,2]
-cbar3 = plt.colorbar(
- a11, ax=ax[1, 0], fraction=0.046, pad=0.04, ticks=ticks_B
-)
+# )
+# cbar5.set_label(r"$\mathrm{log}_{10} \beta \;$")
+ticks_B = [-2, -1, 0, 1, 2]
+cbar3 = plt.colorbar(a11, ax=ax[1, 0], fraction=0.046, pad=0.04, ticks=ticks_B)
cbar3.set_label(r"$\mathrm{log}_{10} B \; [\mu G]$")
-cbar4 = plt.colorbar(
- a11, ax=ax[1, 1], fraction=0.046, pad=0.04, ticks=ticks_B
-)
+cbar4 = plt.colorbar(a11, ax=ax[1, 1], fraction=0.046, pad=0.04, ticks=ticks_B)
cbar4.set_label(r"$\mathrm{log}_{10} B \; [\mu G]$")
ax[2, 1].set_ylabel(r"$y$ [kPc]")
@@ -316,8 +319,8 @@ cbar5.set_label(r"$\mathrm{log}_{10} \frac{h \: \nabla \cdot B}{|B|}$")
cbar6 = plt.colorbar(a21, ax=ax[2, 1], fraction=0.046, pad=0.04)
cbar6.set_label(r"$\mathrm{log}_{10} \frac{h \: \nabla \cdot B}{|B|}$")
-ticksBeta=[0,1]
-cbar7 = plt.colorbar(a30, ax=ax[3, 0], fraction=0.046, pad=0.04,ticks=ticksBeta)
+ticksBeta = [0, 1]
+cbar7 = plt.colorbar(a30, ax=ax[3, 0], fraction=0.046, pad=0.04, ticks=ticksBeta)
cbar7.set_label(r"$\mathrm{log}_{10} \beta $")
# add streamlines
@@ -329,31 +332,31 @@ mass_weighted_vy_map_xy = slice_gas(**common_arguments_xy, project="mass_weighte
vx_map_xy = mass_weighted_vx_map_xy / mass_map_xy
vy_map_xy = mass_weighted_vy_map_xy / mass_map_xy
-vx_map_xy=vx_map_xy.value
-vy_map_xy=vy_map_xy.value
-vplanenorm_map_xy = np.sqrt(vx_map_xy**2+vy_map_xy**2)
-vx_map_xy/=vplanenorm_map_xy
-vy_map_xy/=vplanenorm_map_xy
+vx_map_xy = vx_map_xy.value
+vy_map_xy = vy_map_xy.value
+vplanenorm_map_xy = np.sqrt(vx_map_xy ** 2 + vy_map_xy ** 2)
+vx_map_xy /= vplanenorm_map_xy
+vy_map_xy /= vplanenorm_map_xy
dimx = 512
dimy = 512
-new_x = np.linspace(0,dimx,dimx)
-new_y = np.linspace(0,dimy,dimy)
+new_x = np.linspace(0, dimx, dimx)
+new_y = np.linspace(0, dimy, dimy)
step = 20
-ax[0,1].quiver(
+ax[0, 1].quiver(
new_x[::step],
new_y[::step],
- np.transpose(vx_map_xy.reshape((dimx, dimy)))[::step,::step],
- np.transpose(vy_map_xy.reshape((dimx, dimy)))[::step,::step],
+ np.transpose(vx_map_xy.reshape((dimx, dimy)))[::step, ::step],
+ np.transpose(vy_map_xy.reshape((dimx, dimy)))[::step, ::step],
color="black",
- scale=1.5/step,
- scale_units="xy", # Fixes the arrow length in data coordinates
+ scale=1.5 / step,
+ scale_units="xy", # Fixes the arrow length in data coordinates
pivot="middle"
- #density=1.0,
- #linewidth=0.2,
- #arrowsize=0.4,
- )
+ # density=1.0,
+ # linewidth=0.2,
+ # arrowsize=0.4,
+)
data.gas.mass_weighted_vx = data.gas.masses * v[:, 0]
@@ -363,32 +366,32 @@ mass_weighted_vz_map = slice_gas(**common_arguments, project="mass_weighted_vz")
vx_map = mass_weighted_vx_map / mass_map
vz_map = mass_weighted_vz_map / mass_map
-vx_map=vx_map.value
-vz_map=vz_map.value
-vplanenorm_map = np.sqrt(vx_map**2+vz_map**2)
-vx_map/=vplanenorm_map
-vz_map/=vplanenorm_map
+vx_map = vx_map.value
+vz_map = vz_map.value
+vplanenorm_map = np.sqrt(vx_map ** 2 + vz_map ** 2)
+vx_map /= vplanenorm_map
+vz_map /= vplanenorm_map
dimx = 512
dimy = 512
-new_x = np.linspace(0,dimx,dimx)
-new_y = np.linspace(0,dimy,dimy)
+new_x = np.linspace(0, dimx, dimx)
+new_y = np.linspace(0, dimy, dimy)
step = 20
-ax[0,0].quiver(
+ax[0, 0].quiver(
new_x[::step],
new_y[::step],
- np.transpose(vx_map.reshape((dimx, dimy)))[::step,::step],
- np.transpose(vz_map.reshape((dimx, dimy)))[::step,::step],
+ np.transpose(vx_map.reshape((dimx, dimy)))[::step, ::step],
+ np.transpose(vz_map.reshape((dimx, dimy)))[::step, ::step],
color="black",
- scale=1.5/step,
- scale_units="xy", # Fixes the arrow length in data coordinates
+ scale=1.5 / step,
+ scale_units="xy", # Fixes the arrow length in data coordinates
pivot="middle"
- #density=1.0,
- #linewidth=0.2,
- #arrowsize=0.4,
- )
+ # density=1.0,
+ # linewidth=0.2,
+ # arrowsize=0.4,
+)
data.gas.mass_weighted_Bx = data.gas.masses * B[:, 0]
data.gas.mass_weighted_By = data.gas.masses * B[:, 1]
@@ -396,34 +399,34 @@ mass_weighted_Bx_map_xy = slice_gas(**common_arguments_xy, project="mass_weighte
mass_weighted_By_map_xy = slice_gas(**common_arguments_xy, project="mass_weighted_By")
Bx_map_xy = mass_weighted_Bx_map_xy / mass_map_xy
By_map_xy = mass_weighted_By_map_xy / mass_map_xy
-#Bx_map_xy.convert_to_units(1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s)).values
-#By_map_xy.convert_to_units(1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s)).values
+# Bx_map_xy.convert_to_units(1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s)).values
+# By_map_xy.convert_to_units(1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s)).values
-Bx_map_xy=Bx_map_xy.value
-By_map_xy=By_map_xy.value
-Bplanenorm_map_xy = np.sqrt(Bx_map_xy**2+By_map_xy**2)
-Bx_map_xy/=Bplanenorm_map_xy
-By_map_xy/=Bplanenorm_map_xy
+Bx_map_xy = Bx_map_xy.value
+By_map_xy = By_map_xy.value
+Bplanenorm_map_xy = np.sqrt(Bx_map_xy ** 2 + By_map_xy ** 2)
+Bx_map_xy /= Bplanenorm_map_xy
+By_map_xy /= Bplanenorm_map_xy
dimx = 512
dimy = 512
-new_x = np.linspace(0,dimx,dimx)
-new_y = np.linspace(0,dimy,dimy)
+new_x = np.linspace(0, dimx, dimx)
+new_y = np.linspace(0, dimy, dimy)
step = 20
-ax[1,1].quiver(
+ax[1, 1].quiver(
new_x[::step],
new_y[::step],
- np.transpose(Bx_map_xy.reshape((dimx, dimy)))[::step,::step],
- np.transpose(By_map_xy.reshape((dimx, dimy)))[::step,::step],
+ np.transpose(Bx_map_xy.reshape((dimx, dimy)))[::step, ::step],
+ np.transpose(By_map_xy.reshape((dimx, dimy)))[::step, ::step],
color="black",
- scale=1.5/step,
- scale_units="xy", # Fixes the arrow length in data coordinates
+ scale=1.5 / step,
+ scale_units="xy", # Fixes the arrow length in data coordinates
pivot="middle"
- #density=1.0,
- #linewidth=0.2,
- #arrowsize=0.4,
- )
+ # density=1.0,
+ # linewidth=0.2,
+ # arrowsize=0.4,
+)
data.gas.mass_weighted_Bx = data.gas.masses * B[:, 0]
data.gas.mass_weighted_Bz = data.gas.masses * B[:, 2]
@@ -432,50 +435,52 @@ mass_weighted_Bz_map = slice_gas(**common_arguments, project="mass_weighted_Bz")
Bx_map = mass_weighted_Bx_map / mass_map
Bz_map = mass_weighted_Bz_map / mass_map
-Bx_map=Bx_map.value
-Bz_map=Bz_map.value
-Bplanenorm_map = np.sqrt(Bx_map**2+Bz_map**2)
-Bx_map/=Bplanenorm_map
-Bz_map/=Bplanenorm_map
+Bx_map = Bx_map.value
+Bz_map = Bz_map.value
+Bplanenorm_map = np.sqrt(Bx_map ** 2 + Bz_map ** 2)
+Bx_map /= Bplanenorm_map
+Bz_map /= Bplanenorm_map
dimx = 512
dimy = 512
-new_x = np.linspace(0,dimx,dimx)
-new_y = np.linspace(0,dimy,dimy)
+new_x = np.linspace(0, dimx, dimx)
+new_y = np.linspace(0, dimy, dimy)
step = 20
-ax[1,0].quiver(
+ax[1, 0].quiver(
new_x[::step],
new_y[::step],
- np.transpose(Bx_map.reshape((dimx, dimy)))[::step,::step],
- np.transpose(Bz_map.reshape((dimx, dimy)))[::step,::step],
+ np.transpose(Bx_map.reshape((dimx, dimy)))[::step, ::step],
+ np.transpose(Bz_map.reshape((dimx, dimy)))[::step, ::step],
color="black",
- scale=1.5/step,
- scale_units="xy", # Fixes the arrow length in data coordinates
+ scale=1.5 / step,
+ scale_units="xy", # Fixes the arrow length in data coordinates
pivot="middle"
- #density=1.0,
- #linewidth=0.2,
- #arrowsize=0.4,
- )
+ # density=1.0,
+ # linewidth=0.2,
+ # arrowsize=0.4,
+)
index = np.argmax(normv)
-vpart=normv[index]
+vpart = normv[index]
rpart = r[index]
-Bpart = normB[index]
+Bpart = normB[index]
Berr = np.maximum(h * abs(divB) / normB, 1e-6)
Berrpart = Berr[index]
hpart = h[index]
vpart.convert_to_units(unyt.km / unyt.s)
-rpart.convert_to_units(3.08e18*1e3*unyt.cm)
-hpart.convert_to_units(3.08e18*1e3*unyt.cm)
-Bpart.convert_to_units(1e-7*unyt.g / (unyt.statA * unyt.s * unyt.s))
-print(f'Particle {index}, \nwith maximal velocity {vpart.value} km/s, \nmagnetic field {Bpart.value} muG, \nwith error level {Berrpart.value}, \nwith smoothing length {hpart.value} kPc, \nlocated at {rpart.value} kPc')
+rpart.convert_to_units(3.08e18 * 1e3 * unyt.cm)
+hpart.convert_to_units(3.08e18 * 1e3 * unyt.cm)
+Bpart.convert_to_units(1e-7 * unyt.g / (unyt.statA * unyt.s * unyt.s))
+print(
+ f"Particle {index}, \nwith maximal velocity {vpart.value} km/s, \nmagnetic field {Bpart.value} muG, \nwith error level {Berrpart.value}, \nwith smoothing length {hpart.value} kPc, \nlocated at {rpart.value} kPc"
+)
-part_pixel_coords = 512/2*(1+rpart.value/Lslice_kPc)
+part_pixel_coords = 512 / 2 * (1 + rpart.value / Lslice_kPc)
-#ax[1,0].scatter(part_pixel_coords[0],part_pixel_coords[2],color='red',marker='x')
+# ax[1,0].scatter(part_pixel_coords[0],part_pixel_coords[2],color='red',marker='x')
# add panel with infromation about the run
@@ -510,12 +515,9 @@ ax[3, 1].text(
ax[3, 1].text(
0.5, 0.1, "Physical resistivity $\eta$: $%.2f$ " % (eta), **text_common_args
)
-ax[3, 1].text(
- 0.5, 0.0, "Number of particles $N_p$: $%.0f$ " % (Np), **text_common_args
-)
+ax[3, 1].text(0.5, 0.0, "Number of particles $N_p$: $%.0f$ " % (Np), **text_common_args)
ax[3, 1].axis("off")
fig.tight_layout()
-plt.savefig(sys.argv[2],dpi=220)
-
+plt.savefig(sys.argv[2], dpi=220)
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/run_MWsize.sh b/examples/MHDTests/CoolingHaloWithSpin_MHD/run.sh
similarity index 58%
rename from examples/MHDTests/CoolingHaloWithSpin_MHD/run_MWsize.sh
rename to examples/MHDTests/CoolingHaloWithSpin_MHD/run.sh
index 26c5263416dac7d3e02ecde26f81c8e2fb017cac..cfbfb42dd8e7aa3e123e13aa97f276c129d16802 100755
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/run_MWsize.sh
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/run.sh
@@ -2,7 +2,7 @@
# Generate the initial conditions if they are not present.
echo "Generating initial conditions for the MHD isothermal potential box example..."
-python3 makeIC_MWsize.py 1000000
+python3 makeIC_MWsize.py 45000
# Run SWIFT with external potential, SPH and cooling
-../../../swift --external-gravity --hydro --cooling --threads=70 cooling_halo_MWsize.yml 2>&1 | tee output.log
+../../../swift --external-gravity --self-gravity --hydro --cooling --threads=4 cooling_halo.yml 2>&1 | tee output.log
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/run_dwarfsize.sh b/examples/MHDTests/CoolingHaloWithSpin_MHD/run_dwarfsize.sh
deleted file mode 100755
index 32f5535e6b30ff94f66cd29711b6f2228b7cbbae..0000000000000000000000000000000000000000
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/run_dwarfsize.sh
+++ /dev/null
@@ -1,8 +0,0 @@
-#!/bin/bash
-
-# Generate the initial conditions if they are not present.
-echo "Generating initial conditions for the MHD isothermal potential box example..."
-python3 makeIC_dwarfsize.py 1000000
-
-# Run SWIFT with external potential, SPH and cooling
-../../../swift --external-gravity --hydro --cooling --threads=70 cooling_halo_dwarfsize.yml 2>&1 | tee output.log
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/statistics_reader.py b/examples/MHDTests/CoolingHaloWithSpin_MHD/statistics_reader.py
index 42ef23274d2d859b6a9155a71f06e885b371487f..63dd09c3d4ba6d727c5c94b6e169283ad5c58122 100644
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/statistics_reader.py
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/statistics_reader.py
@@ -7,6 +7,7 @@ import sys
from glob import glob
from swiftsimio import load_statistics
import unyt as u
+
filename = sys.argv[1]
outputname = sys.argv[2]
@@ -19,18 +20,18 @@ Eint = data.int_energy.to(u.erg)
Emag = data.mag_energy.to(u.erg)
Epot = np.abs(data.pot_energy.to(u.erg))
-Etot = np.abs(Ekin+Eint+Emag-Epot)
+Etot = np.abs(Ekin + Eint + Emag - Epot)
-name = ''
+name = ""
-#print(max(Epot))
+# print(max(Epot))
fig, ax = plt.subplots(1, 1, sharex=True, figsize=(7, 5))
-ax.plot(time, Emag,label='Emag'+name,color='green',linestyle='dashed')
-ax.plot(time, Ekin,label='Ekin'+name,color='blue',linestyle='dashed')
-ax.plot(time, Eint,label='Eint'+name, color='red',linestyle='dashed')
-ax.plot(time, Epot,label='|Epot|'+name,color='brown', linestyle='dotted')
-#ax.plot(time, Epot,label='Etot'+name,color='black', linestyle='dotted')
+ax.plot(time, Emag, label="Emag" + name, color="green", linestyle="dashed")
+ax.plot(time, Ekin, label="Ekin" + name, color="blue", linestyle="dashed")
+ax.plot(time, Eint, label="Eint" + name, color="red", linestyle="dashed")
+ax.plot(time, Epot, label="|Epot|" + name, color="brown", linestyle="dotted")
+# ax.plot(time, Epot,label='Etot'+name,color='black', linestyle='dotted')
ax.set_xlabel("time [Gyr]")
ax.set_ylabel(f"Energy [erg]")
ax.set_ylim([1e48, 1e60])
@@ -38,4 +39,3 @@ ax.set_yscale("log")
ax.grid()
ax.legend()
plt.savefig(outputname, dpi=220)
-