diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/1Dprofiles.py b/examples/MHDTests/CoolingHaloWithSpin_MHD/1Dprofiles.py
index ee568c98ee39df18e615f2c750589a68f2f81294..f47a488a905f39ce0e38149d1d5e30241c618f70 100644
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/1Dprofiles.py
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/1Dprofiles.py
@@ -18,8 +18,8 @@ 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_cgs = 100.0 * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
+h_dim = 0.704 #0.681 # 0.67777 # hubble parameter
+H_0_cgs = 100.0 * h_dim * 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
@@ -28,7 +28,7 @@ H_0_cgs = 100.0 * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
f_b = 0.17
c_200 = 7.2
spin_lambda = 0.05
-nH_max_cgs = 1e0
+nH_max_cgs = 1e2
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)
@@ -79,76 +79,43 @@ n_gas = data.metadata.n_gas
rho = data.gas.densities
m = data.gas.masses
coords = data.gas.coordinates
+coords_center = coords - data.metadata.boxsize / 2
+radius = np.linalg.norm(coords_center,axis=1)
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))
+
+# calculate specific angular momentum j(r)
+rotation_axis = np.array([0,0,1])
+e_phi = np.cross(rotation_axis[None,:],coords_center)
+e_phi = e_phi/np.linalg.norm(e_phi,axis=1)[:,None]
+v_phi = velocities[:,0]*e_phi[:,0]+velocities[:,1]*e_phi[:,1]+velocities[:,2]*e_phi[:,2]
+e_r = np.cross(e_phi,rotation_axis[None,:])
+axis_dist = coords_center[:,0]*e_r[:,0]+coords_center[:,1]*e_r[:,1]+coords_center[:,2]*e_r[:,2]
+j = v_phi * radius**2 / axis_dist
+omega = v_phi / axis_dist
+Jtot = np.sum(omega * axis_dist**2 * m)
P = data.gas.pressures
T = data.gas.temperatures
-
B = data.gas.magnetic_flux_densities
-
Bx, By = B[:, 0], B[:, 1]
-
normB = np.sqrt(B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2)
-
divB = data.gas.magnetic_divergences
-
h = data.gas.smoothing_lengths
-
errB = np.log10(h * abs(divB) / normB)
-# Generate mass weighted maps of quantities of interest
-data.gas.mass_weighted_densities = data.gas.masses * rho
-
-data.gas.mass_weighted_pressures = data.gas.masses * P
-
-data.gas.mass_weighted_Bx = data.gas.masses * Bx
-data.gas.mass_weighted_By = data.gas.masses * By
+rho_units = unyt.g * unyt.cm ** (-3)
+r_units = 1e3 * PARSEC_IN_CGS * unyt.cm
+P_units = MSOL_IN_CGS / PARSEC_IN_CGS / GYR_IN_CGS ** 2 * (unyt.g / (unyt.cm * unyt.s**2))
+j_units = (1e3 * PARSEC_IN_CGS) / GYR_IN_CGS * (unyt.cm**2/unyt.s)
+nH_units = unyt.cm ** (-3)
-data.gas.mass_weighted_errB = data.gas.masses * errB
-
-data.gas.mass_weighted_T = data.gas.masses * T
-
-common_arguments = dict(
- data=data, z_slice=0.5 * data.metadata.boxsize[2], resolution=512, parallel=True
-)
-
-mass_map = slice_gas(**common_arguments, project="masses")
-
-mass_weighted_density_map = slice_gas(
- **common_arguments, project="mass_weighted_densities"
-)
-
-mass_weighted_pressure_map = slice_gas(
- **common_arguments, project="mass_weighted_pressures"
-)
-
-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_errB_map = slice_gas(**common_arguments, project="mass_weighted_errB")
-mass_weighted_T_map = slice_gas(**common_arguments, project="mass_weighted_T")
-
-# Take out mass dependence
-density_map = mass_weighted_density_map / mass_map
-pressure_map = mass_weighted_pressure_map / mass_map
-Bx_map = mass_weighted_Bx_map / mass_map
-By_map = mass_weighted_By_map / mass_map
-errB_map = mass_weighted_errB_map / mass_map
-T_map = mass_weighted_T_map / mass_map
-
-map_pixel_length = len(mass_map)
-
-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))
+rho.convert_to_units(rho_units)
+P.convert_to_units(P_units)
+j.convert_to_units(j_units)
+n_H = rho/(m_H_cgs * unyt.g)
+radius = np.linalg.norm(coords_center,axis=1)
# Plot maps
plt.rcParams.update({"font.size": 16})
@@ -158,42 +125,19 @@ 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)
-slice_ind = int(np.floor(y0 * map_pixel_length))
-
# plot density
-axs[0].plot(x, n_H[:, slice_ind], "k-", color="black")
+axs[0].scatter(radius.to(r_units), n_H.to(nH_units), s=0.1, 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")
+# plot Pressure
+axs[1].scatter(radius.to(r_units), P.to(P_units), s=0.1, 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)
-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)
-
# 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)
-
+axs[2].scatter(radius.to(r_units), j.to(j_units), s=0.1, color='black')
+axs[2].set_yscale("log")
# NFW-like gas density profile
def rho_r(r_value, f_b, M_200_cgs, r_200_cgs, c_200):
@@ -273,100 +217,116 @@ def T_vs_r(P0_cgs, r_value, r_max, f_b, M_200_cgs, r_200_cgs, c_200):
)
return result_cgs
+# 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 jsimple(r_value, j_max, s, f_b, M_200_cgs, r_200_cgs, c_200):
+ return (
+ j_max
+ * r_value**s
+ )
-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
-)
+rmax = (r_200_cgs*unyt.cm).to(radius.units) #np.max(radius)
+rmin = np.min(radius)
+r_analytic = np.logspace(np.log10(rmin.value),np.log10(rmax.value),1000)*radius.units
+#r_analytic = np.linspace(rmin.value,rmax.value,10000)*radius.units
+rho_analytic = rho_r((r_analytic/(r_200_cgs*unyt.cm)).value, f_b, M_200_cgs, r_200_cgs, c_200)
+Mgas_analytic = Mgas_r((r_analytic/(r_200_cgs*unyt.cm)).value, f_b, M_200_cgs, r_200_cgs, c_200)
+nH_analytic = rho_analytic/m_H_cgs * nH_units
-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
+ (rmax/(r_200_cgs*unyt.cm)).value, f_b, M_200_cgs, r_200_cgs, c_200
+)
-Pressure_UNIT = MSOL_IN_CGS / PARSEC_IN_CGS / GYR_IN_CGS ** 2
-P_analytic = [
+P0_cgs = rho0_cgs * kb_cgs * T0_cgs / m_H_cgs # gas pressure on the edge of the box
+P_analytic =np.array( [
P_vs_r(
P0_cgs,
- np.abs(x[i] / r_200_kpc),
- Lbox_kpc / r_200_kpc * np.sqrt(3) / 2,
+ (r_analytic[i]/(r_200_cgs*unyt.cm)).value,
+ (rmax/(r_200_cgs*unyt.cm)).value,
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
+ for i in range(len(r_analytic))
+]) * (unyt.g/(unyt.cm*unyt.s**2)) # gas particle internal energies
+# Normalize to Bullock
+mask_r200 = radius<=(r_200_cgs * unyt.cm)
+Jtot = np.sum(omega[mask_r200] * axis_dist[mask_r200]**2 * m[mask_r200])
+jsp = Jtot / np.sum(m[mask_r200])
-# 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)
+jsp_B = (np.sqrt(2) * v_200_cgs * (unyt.cm/unyt.s))*r_200_cgs*(unyt.cm)*spin_lambda
+print('Total spin ratio: ',jsp_B.to(j_units).value/jsp.to(j_units).value)
-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_ylabel(r"$n_H(r)$ $[cm^{-3}]$")
-axs[0].axvline(x=r_200_kpc, color="gray", ls="dashed", label="$R_200$")
+############
+
+j_1_analytic = j((r_analytic/(r_200_cgs*unyt.cm)).value, 1, 1, f_b, M_200_cgs, r_200_cgs, c_200)
+j_1_analytic_simple = jsimple((r_analytic/(r_200_cgs*unyt.cm)).value, 1, 2, f_b, M_200_cgs, r_200_cgs, c_200)
+
+int_dOmega = 8*np.pi/3 # integral over angles from axial distance squared: int Sin^2 (theta)dphi d Cos (theta)
+
+jsp_1_analytic = (np.sum(j_1_analytic[:-1] * rho_analytic[:-1] * (unyt.g/unyt.cm**3) * int_dOmega * r_analytic[:-1]**2*np.diff(r_analytic)))/(np.sum(rho_analytic[:-1] * (unyt.g/unyt.cm**3) * 4 * np.pi * r_analytic[:-1]**2*np.diff(r_analytic)))
+jsp_1_analytic_simple = (np.sum(j_1_analytic_simple[:-1] * rho_analytic[:-1] * (unyt.g/unyt.cm**3) * int_dOmega * r_analytic[:-1]**2*np.diff(r_analytic)))/(np.sum(rho_analytic[:-1] * (unyt.g/unyt.cm**3) * 4 * np.pi * r_analytic[:-1]**2*np.diff(r_analytic)))
+
+j_analytic = j_1_analytic / jsp_1_analytic * jsp_B
+
+j_analytic_simple = j_1_analytic_simple / jsp_1_analytic_simple * jsp_B
+
+
+
+#
+# plot density
+axs[0].plot(r_analytic.to(r_units), nH_analytic.to(nH_units), color='red',label = r'$\rho_{\rm NFW}^{gas}(r)$')
+axs[0].set_yscale("log")
axs[0].set_xscale("log")
-axs[0].set_xlim(0.5, 300)
+axs[0].set_xlim([rmin,rmax])
+axs[0].legend()
-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_ylabel(r"$P(r)$ $[M_{sol} \cdot pc^{-1} \cdot Gyr^{-2}]$")
+# plot Pressure
+axs[1].plot(r_analytic.to(r_units), P_analytic.to(P_units), color='red',label = r'$P_{hyd. eq.}(r)$')
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].set_xlim([rmin,rmax])
+axs[1].legend()
-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_ylabel(r"$M_{gas}(<r)$ $[M_{sol}]$")
-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$"
-)
+# plot j(r)
+axs[2].plot(r_analytic.to(r_units), j_analytic.to(j_units), color='red',label = r'$\rm j(r)\sim M(r)$')
+axs[2].plot(r_analytic.to(r_units), j_analytic_simple.to(j_units), color='red', linestyle='dashed',label=r'$\rm j(r)\sim r^2$')
axs[2].set_xscale("log")
-axs[2].set_xlim(0.5, 300)
+axs[2].set_xlim([rmin,rmax])
+axs[2].legend()
+
+# limits
+axs[0].set_ylim(1e-7, 1e2)
+axs[0].set_yticks(np.logspace(-7, 2, 10))
+axs[1].set_ylim(1e1, 1e9)
+axs[1].set_yticks(np.logspace(1, 9, 9))
+axs[2].set_ylim(1e18, 1e27)
+axs[2].set_yticks(np.logspace(18, 27, 10))
+
+axs[0].set_xlim(1e-2, 1e3)
+axs[1].set_xlim(1e-2, 1e3)
+axs[2].set_xlim(1e-2, 1e3)
+# add labels
+axs[0].set_ylabel(r"$n_H(r)$ $[cm^{-3}]$")
+axs[1].set_ylabel(r"$P(r)$ $[\rm M_{sol} \cdot pc^{-1} \cdot Gyr^{-2}]$")
+axs[2].set_ylabel(r"$j(r)$ $[\rm kpc^{2} \cdot Gyr^{-1} ]$")
+
+# vertical lines
+axs[0].axvline(x=(r_200_cgs*unyt.cm).to(r_units),color='gray',label=r'$R_{200}$')
+axs[1].axvline(x=(r_200_cgs*unyt.cm).to(r_units),color='gray',label=r'$R_{200}$')
+axs[2].axvline(x=(r_200_cgs*unyt.cm).to(r_units),color='gray',label=r'$R_{200}$')
Ninfo = 3
text_common_args = dict(
@@ -376,7 +336,7 @@ text_common_args = dict(
axs[Ninfo].text(
0.5,
0.8,
- r"Cooling Halo with Spin at $t=%.2f$, $y_0 = %.4f$" % (data.metadata.time, y0),
+ r"Cooling Halo with Spin, IC quality check.",
**text_common_args,
)
axs[Ninfo].text(0.5, 0.7, r"SWIFT %s" % git.decode("utf-8"), **text_common_args)
@@ -389,16 +349,7 @@ axs[Ninfo].text(
**text_common_args,
)
axs[Ninfo].text(
- 0.5, 0.3, r"Artificial diffusion: $%.2f$ " % (artDiffusion), **text_common_args
-)
-axs[Ninfo].text(
- 0.5,
- 0.2,
- r"Dedner Hyp, Hyp_div(v), Par: $%.2f,%.2f,%.2f$ " % (dedHyp, dedHypDivv, dedPar),
- **text_common_args,
-)
-axs[Ninfo].text(
- 0.5, 0.1, r"Physical resistivity: $\eta$: $%.2f$ " % (eta), **text_common_args
+ 0.5, 0.3, r"Dimenstionless hubble constant: $%.3f$ " % h_dim, **text_common_args
)
axs[Ninfo].text(
0.5, 0.0, r"Number of particles: $N_p$: $%.0f$ " % (n_gas), **text_common_args
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/README b/examples/MHDTests/CoolingHaloWithSpin_MHD/README
index adc8ad69e6a50d981c11eb868785b30a1be9e552..6d246217c5ef7e6d0b6343ccb797bbfd5824866b 100644
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/README
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/README
@@ -24,15 +24,22 @@ collapse into a spinning disc.
Compilation
-----------
-for runs with EAGLE cooling+chemistry+entropy floor
-./configure --with-spmhd=direct-induction --with-cooling=EAGLE --with-chemistry=EAGLE --with-entropy-floor=EAGLE --with-ext-potential=nfw --with-kernel=quintic-spline --disable-hand-vec --disable-compiler-warnings --disable-doxygen-doc
+
+for runs with EAGLE cooling+chemistry+entropy floor+star formation + feedback
+./configure --with-spmhd=direct-induction --with-cooling=EAGLE --with-chemistry=EAGLE --with-entropy-floor=EAGLE --with-ext-potential=nfw --with-kernel=quintic-spline --disable-hand-vec --disable-compiler-warnings --disable-doxygen-doc --with-stars=EAGLE --with-star-formation=EAGLE --with-tracers=EAGLE --with-feedback=EAGLE
+
+for runs with EAGLE cooling+chemistry+entropy floor+star formation
+./configure --with-spmhd=direct-induction --with-cooling=EAGLE --with-chemistry=EAGLE --with-entropy-floor=EAGLE --with-ext-potential=nfw --with-kernel=quintic-spline --disable-hand-vec --disable-compiler-warnings --disable-doxygen-doc --with-stars=EAGLE --with-star-formation=EAGLE
for runs with lambda cooling
./configure --with-spmhd=direct-induction --with-cooling=const-lambda --with-ext-potential=nfw --with-kernel=quintic-spline --disable-hand-vec --disable-compiler-warnings --disable-doxygen-doc
+Note: when running with self gravity but without feedback excessive clumping will happen. Since there is no density limiter, the clumps collapse until they reach sizes ~ softening length. Rotation can dilute and destabilize these clums, so adding more angular momentum might help in this case.
+
Setup
-----
-We follow R. Pakmor, V. Springel, 1212.1452
+IC generation similar to R. Pakmor, V. Springel, 1212.1452
+Run parameters for EAGLE model mostly follow run L025N0376 from J. Schaye et al. 1407.7040
Snapshot slices
---------------
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/cooling_halo.yml b/examples/MHDTests/CoolingHaloWithSpin_MHD/cooling_halo.yml
index 6b8b89af1cc7cd641b6a819a8ca4eaae612b5db3..c875f3677dc51720e053cc543b74e594799f86e7 100644
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/cooling_halo.yml
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/cooling_halo.yml
@@ -1,80 +1,126 @@
-# Define the system of units to use internally.
+# 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
+ UnitMass_in_cgs: 1.9891E43 # 10^10 solar masses
+ UnitLength_in_cgs: 3.08567758E21 # 1 kpc
+ UnitVelocity_in_cgs: 1E5 # km/s
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).
+ time_end: 3. # 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
+ delta_time: 0.01023 # Time difference between consecutive outputs (in internal units) 0.001023 TU = 10 Myr
+ compression: 4 # Compress the snapshots
+ recording_triggers_part: [-1, -1] # Not recording as we have many snapshots
+ recording_triggers_bpart: [-1, -1] # Not recording as we have many snapshots
# 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)
-
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 0.01023 # Time difference between consecutive outputs (in internal units)
+ compression: 4 # Compress the snapshots
+ recording_triggers_part: [-1, -1] # Not recording as we have many snapshots
+ recording_triggers_bpart: [-1, -1] # Not recording as we have many snapshots
+
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.595 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
CFL_condition: 0.05 # Courant-Friedrich-Levy condition for time integration.
minimal_temperature: 1e2 # Kelvin
- h_tolerance: 1e-6
-
+ h_tolerance: 1e-6 # smoothing length accuracy
+ h_max: 10.
+
# Parameters for MHD schemes
MHD:
monopole_subtraction: 1.0
artificial_diffusion: 0.0
- hyperbolic_dedner: 1.0 #1.0
- hyperbolic_dedner_divv: 0.0
+ hyperbolic_dedner: 1.0
+ hyperbolic_dedner_divv: 0.5
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
-
-# Note: softening 2.6 below is for 1.4e7 Msol particle mass to impose density cut 1e2 cm^-3
+ periodic: 0
+ stars_smoothing_length: 0.5
# 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: 1000000 # M200 of the galaxy disk
- h: 0.681 # reduced Hubble constant (using FLAMINGO cosmology) (value does not specify the used units!)
+ M_200: 100 # M200 of the galaxy disk
+ h: 0.704 # reduced Hubble constant (using FLAMINGO cosmology) (value does not specify the used units!)
concentration: 7.2 # 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.569 #2.6 # Softening size (internal units)
- #diskfraction: 0.0 # Disk mass fraction
- #bulgefraction: 0.0 # Bulge mass fraction
+ epsilon: 0.7 # Softening size (internal units)
# 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.569 #2.6 # Physical softening length (in internal units). Depends on resolution!
- mesh_side_length: 64
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ MAC: adaptive
+ theta_cr: 0.6 # Opening angle (Multipole acceptance criterion).
+ epsilon_fmm: 0.001
+ max_physical_DM_softening: 0.7 # Physical softening length (in internal units).
+ max_physical_baryon_softening: 0.7 # Physical softening length (in internal units).
+
+# Parameters for the stars
+Stars:
+ overwrite_birth_time: 1 # Make sure the stars in the ICs do not do any feedback
+ birth_time: -1. # by setting all of their birth times to -1
+ timestep_age_threshold_Myr: 10. # Age at which stars switch from young to old (in Mega-years).
+ max_timestep_young_Myr: 0.1 # Maximal time-step length of young stars (in Mega-years).
+ max_timestep_old_Myr: 1.0 # Maximal time-step length of old stars (in Mega-years).
+ luminosity_filename: ./photometry
+
+# Standard EAGLE cooling options
+EAGLECooling:
+ dir_name: ./coolingtables/ # Location of the Wiersma+09 cooling tables
+ H_reion_z: 11.5 # Redshift of Hydrogen re-ionization
+ H_reion_eV_p_H: 2.0
+ 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
-# 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
+# Solar abundances
+EAGLEChemistry:
+ init_abundance_metal: 0.0129 # Inital fraction of particle mass in *all* metals
+ init_abundance_Hydrogen: 0.7065 # Inital fraction of particle mass in Hydrogen
+ init_abundance_Helium: 0.2806 # Inital fraction of particle mass in Helium
+ init_abundance_Carbon: 0.00207 # Inital fraction of particle mass in Carbon
+ init_abundance_Nitrogen: 0.000836 # Inital fraction of particle mass in Nitrogen
+ init_abundance_Oxygen: 0.00549 # Inital fraction of particle mass in Oxygen
+ init_abundance_Neon: 0.00141 # Inital fraction of particle mass in Neon
+ init_abundance_Magnesium: 0.000591 # Inital fraction of particle mass in Magnesium
+ init_abundance_Silicon: 0.000683 # Inital fraction of particle mass in Silicon
+ init_abundance_Iron: 0.0011 # Inital fraction of particle mass in Iron
+# EAGLE star formation parameters
+EAGLEStarFormation:
+ SF_threshold: Zdep # Zdep (Schaye 2004) or Subgrid
+ SF_model: PressureLaw # PressureLaw (Schaye et al. 2008) or SchmidtLaw
+ gas_fraction: 0.3 # The gas fraction used internally by the model.
+ KS_normalisation: 1.515e-4 # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+ KS_exponent: 1.4 # The exponent of the Kennicutt-Schmidt law.
+ min_over_density: 100.0 # The over-density above which star-formation is allowed.
+ KS_high_density_threshold_H_p_cm3: 1e3 # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+ KS_high_density_exponent: 2.0 # Slope of the Kennicut-Schmidt law above the high-density threshold.
+ EOS_entropy_margin_dex: 0.3 # When using Z-based SF threshold, logarithm base 10 of the maximal entropy above the EOS at which stars can form.
+ threshold_norm_H_p_cm3: 0.1 # When using Z-based SF threshold, normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+ threshold_Z0: 0.002 # When using Z-based SF threshold, reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+ threshold_slope: -0.64 # When using Z-based SF threshold, slope of the metal-dependant star formation threshold
+ threshold_max_density_H_p_cm3: 10.0 # When using Z-based SF threshold, maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+ threshold_temperature1_K: 1000 # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming.
+ threshold_temperature2_K: 31622 # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming if also above the density limit.
+ threshold_number_density_H_p_cm3: 10 # When using subgrid-based SF threshold, subgrid number density above which gas is star-forming if also below the second temperature limit.
+ #
# 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.
@@ -86,39 +132,46 @@ EAGLEEntropyFloor:
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 primordial abundances
-EAGLEChemistry:
- init_abundance_metal: 0.
- init_abundance_Hydrogen: 0.755
- init_abundance_Helium: 0.245
- init_abundance_Carbon: 0.
- init_abundance_Nitrogen: 0.
- init_abundance_Oxygen: 0.
- init_abundance_Neon: 0.
- init_abundance_Magnesium: 0.
- init_abundance_Silicon: 0.
- init_abundance_Iron: 0.
-
-# 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
-
+# EAGLE feedback model with constant feedback energy fraction
+EAGLEFeedback:
+ use_SNII_feedback: 1 # Global switch for SNII thermal (stochastic) feedback.
+ use_SNIa_feedback: 1 # Global switch for SNIa thermal (continuous) feedback.
+ use_AGB_enrichment: 1 # Global switch for enrichement from AGB stars.
+ use_SNII_enrichment: 1 # Global switch for enrichement from SNII stars.
+ use_SNIa_enrichment: 1 # Global switch for enrichement from SNIa stars.
+ filename: ./yieldtables/ # Path to the directory containing the EAGLE yield tables.
+ IMF_min_mass_Msun: 0.1 # Minimal stellar mass considered for the Chabrier IMF in solar masses.
+ IMF_max_mass_Msun: 100.0 # Maximal stellar mass considered for the Chabrier IMF in solar masses.
+ SNII_min_mass_Msun: 8.0 # Minimal mass considered for SNII stars in solar masses.
+ SNII_max_mass_Msun: 100.0 # Maximal mass considered for SNII stars in solar masses.
+ SNII_feedback_model: MinimumDistance # Feedback modes: Random, Isotropic, MinimumDistance, MinimumDensity
+ SNII_sampled_delay: 1 # Sample the SNII lifetimes to do feedback.
+ SNII_delta_T_K: 3.16228e7 # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
+ SNII_delta_v_km_p_s: 200 # Velocity kick applied by the stars when doing SNII feedback (in km/s).
+ SNII_energy_erg: 1.0e51 # Energy of one SNII explosion in ergs.
+ SNII_energy_fraction_function: EAGLE # Type of functional form to use for scaling the energy fraction with density and metallicity ('EAGLE', 'Separable', or 'Independent').
+ SNII_energy_fraction_min: 1.0 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_max: 1.0 # Maximal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_Z_0: 0.0012663729 # Pivot point for the metallicity dependance of the SNII energy fraction (metal mass fraction).
+ SNII_energy_fraction_n_0_H_p_cm3: 1.4588 # Pivot point for the birth density dependance of the SNII energy fraction in cm^-3.
+ SNII_energy_fraction_n_Z: 0.8686 # Power-law for the metallicity dependance of the SNII energy fraction.
+ SNII_energy_fraction_n_n: 0.8686 # Power-law for the birth density dependance of the SNII energy fraction.
+ SNII_energy_fraction_use_birth_density: 0 # Are we using the density at birth to compute f_E or at feedback time?
+ SNII_energy_fraction_use_birth_metallicity: 0 # Are we using the metallicity at birth to compuote f_E or at feedback time?
+ SNIa_DTD: Exponential # Functional form of the SNIa delay time distribution.
+ SNIa_DTD_delay_Gyr: 0.04 # Stellar age after which SNIa start in Gyr (40 Myr corresponds to stars ~ 8 Msun).
+ SNIa_DTD_exp_timescale_Gyr: 2.0 # Time-scale of the exponential decay of the SNIa rates in Gyr.
+ SNIa_DTD_exp_norm_p_Msun: 0.002 # Normalisation of the SNIa rates in inverse solar masses.
+ SNIa_energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
+ AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
+ stellar_evolution_age_cut_Gyr: 0.1 # Age in Gyr above which the enrichment is downsampled.
+ stellar_evolution_sampling_rate: 10 # Number of time-steps in-between two enrichment events for a star above the age threshold.
+ SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
+ SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
+ SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
+ SNII_yield_factor_Nitrogen: 1.0 # (Optional) Correction factor to apply to the Nitrogen yield from the SNII channel.
+ SNII_yield_factor_Oxygen: 1.0 # (Optional) Correction factor to apply to the Oxygen yield from the SNII channel.
+ SNII_yield_factor_Neon: 1.0 # (Optional) Correction factor to apply to the Neon yield from the SNII channel.
+ SNII_yield_factor_Magnesium: 2.0 # (Optional) Correction factor to apply to the Magnesium yield from the SNII channel.
+ SNII_yield_factor_Silicon: 1.0 # (Optional) Correction factor to apply to the Silicon yield from the SNII channel.
+ SNII_yield_factor_Iron: 0.5 # (Optional) Correction factor to apply to the Iron yield from the SNII channel.
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/getEaglePhotometryTable.sh b/examples/MHDTests/CoolingHaloWithSpin_MHD/getEaglePhotometryTable.sh
new file mode 100755
index 0000000000000000000000000000000000000000..ee9c3b422f19518612416da0913b162fd4a120ff
--- /dev/null
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/getEaglePhotometryTable.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/YieldTables/EAGLE/photometry.tar.gz
+tar -xf photometry.tar.gz
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/getYieldTable.sh b/examples/MHDTests/CoolingHaloWithSpin_MHD/getYieldTable.sh
new file mode 100755
index 0000000000000000000000000000000000000000..26eef020cab82acee2c80e88089df1790b281eab
--- /dev/null
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/getYieldTable.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/YieldTables/EAGLE/yieldtables.tar.gz
+tar -xf yieldtables.tar.gz
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC.py b/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC.py
index ecd957484ca3440b16b2d42c9add49b0025e5d56..a06466e55816062f2c98cc4ca588f99d60adb6ee 100644
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC.py
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC.py
@@ -41,30 +41,36 @@ 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.704 # 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
# 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-3 micro Gauss
B0_cgs = np.sqrt(CONST_MU0_CGS / (4.0 * np.pi)) * B0_Gaussian_Units
# SPH
-eta = 1.3663 # kernel smoothing
+eta = 1.595 # kernel smoothing
# Additional options
spin_lambda_choice = (
- "Bullock"
+ "Peebles"
) # which definition of spin parameter to use (Peebles or Bullock)
save_to_vtk = False
plot_v_distribution = False
+AMP = "r^s" # sets angular momentum profile to be a function of M(r) or r^s
+AMPs = 2 # power in agnular momentum profile
+
+with_noise = True # add gaussian uncertainty for particle positions
+noise_sigma = 0.5
+noise_mean = 0.0
+
# From this we can find the virial radius, the radius within which the average density of the halo is
# 200. * the mean matter density
@@ -112,9 +118,11 @@ print("G=", const_G)
# Parameters
periodic = 1 # 1 For periodic box
-boxSize = 4.0
+boxSize = 2.0 # in units of r_200
+R_max = boxSize/6 #boxSize/6 # set maximal halo radius (we simulate only part of a halo)
G = const_G
N = int(sys.argv[1]) # Number of particles
+N = int(N*6/np.pi) # renormalize number of particles to get required N after cutting a sphere
IAsource = "IAfile"
# Create the file
@@ -190,10 +198,9 @@ coords -= 0.5
# cut a sphere
r_uni = np.linalg.norm(coords, axis=1)
coords = coords[r_uni <= 0.5]
-coords *= boxSize * np.sqrt(3)
+coords *= 2*R_max
# calculate max distance from the center (units of r_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)
@@ -301,6 +308,18 @@ r = np.logspace(
round(10 * N ** (1 / 3)),
)
+# Add noise
+if with_noise:
+ radius = np.sqrt(
+ (coords[:, 0]) ** 2
+ + (coords[:, 1]) ** 2
+ + (coords[:, 2]) ** 2
+ )
+ l = (gas_particle_mass*M_200_cgs/rho_r(radius, f_b, M_200_cgs, r_200_cgs, c_200))**(1/3)/r_200_cgs
+ random_amps = noise_sigma * np.random.randn(len(coords),3) + noise_mean
+ coords += random_amps * l[:,None]
+
+
# 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])))
@@ -397,11 +416,17 @@ print("j_sp =", j_sp)
# is proportional to j/r^2
# 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 j(r_value, j_max, s, f_b, M_200_cgs, r_200_cgs, c_200, angular_momentum_profile="M(r)"):
+ if angular_momentum_profile == "M(r)":
+ return (
+ j_max
+ * (Mgas_r(r_value, f_b, M_200_cgs, r_200_cgs, c_200) / (M_200_cgs * f_b)) ** s
+ )
+ elif angular_momentum_profile=="r^s": # for rigid body rotation use r^2
+ return (
+ j_max
+ * r_value**s
+ )
omega = np.zeros((N, 3))
@@ -410,7 +435,7 @@ radius_vect = 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
+ omega[i, 2] = j(radius[i], 1, AMPs, f_b, M_200_cgs, r_200_cgs, c_200,AMP) / 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, :])
@@ -418,7 +443,8 @@ for i in range(N):
# select particles inside R_200
mask_r_200 = radius <= 1
Lp_tot = np.sum(l[mask_r_200], axis=0)
-Lp_tot_abs = np.linalg.norm(Lp_tot)
+Lp_tot_abs = np.linalg.norm(Lp_tot)
+jsp = Lp_tot_abs/gas_mass
normV = np.linalg.norm(v, axis=1)
if spin_lambda_choice == "Peebles":
@@ -430,17 +456,15 @@ if spin_lambda_choice == "Peebles":
/ (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)
+ jsp_P = (
+ const_G * (gas_mass) ** (3 / 2) / np.sqrt(np.abs(Ep_tot))
)
- v *= spin_lambda / calc_spin_par_P
- l *= spin_lambda / calc_spin_par_P
+ v *= jsp_P/jsp
+
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
+ jsp_B = (np.sqrt(2) * v_200 )*spin_lambda
+ v *= jsp_B/jsp
# Draw velocity distribution
if plot_v_distribution:
@@ -488,9 +512,8 @@ 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
+
+l = (gas_particle_mass*M_200_cgs/rho_r(radius, f_b, M_200_cgs, r_200_cgs, c_200))**(1/3)/r_200_cgs
h = np.full((N,), eta * l)
ds = grp.create_dataset("SmoothingLength", (N,), "f")
ds[()] = h
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/plotSolutionReference.py b/examples/MHDTests/CoolingHaloWithSpin_MHD/plotSolutionReference.py
index a563d432517d48f16a5a76f48af9729ac6138c01..efd156b9cab4dd41b251bf6f27613e9a08c849c2 100644
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/plotSolutionReference.py
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/plotSolutionReference.py
@@ -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
Lslice = Lslice_kPc * 3.08e18 * 1e3 * unyt.cm
visualise_region_xy = [
@@ -204,7 +204,7 @@ fig, ax = plt.subplots(ny, nx, sharey=True, figsize=(5 * nx, 5 * ny))
a00 = ax[0, 0].contourf(
np.log10(nH_map.value),
# np.log10(density_map.value),
- cmap="jet", # "plasma",#"gist_stern",#"RdBu",
+ cmap="plasma",
extend="both",
levels=np.linspace(-4, 2, 100),
# levels=np.linspace(-7, -1, 100),
@@ -213,7 +213,7 @@ a00 = ax[0, 0].contourf(
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",
+ cmap="plasma",
extend="both",
levels=np.linspace(-4, 2, 100),
# levels=np.linspace(-7, -1, 100),
@@ -236,13 +236,13 @@ a01 = ax[0, 1].contourf(
# )
a10 = ax[1, 0].contourf(
np.maximum(np.log10(normB_map.value), -6),
- cmap="jet", # "viridis", #"nipy_spectral_r",
+ cmap="viridis",
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="viridis",
extend="both",
levels=np.linspace(-2, 2, 100),
)
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/plotSpectrum.py b/examples/MHDTests/CoolingHaloWithSpin_MHD/plotSpectrum.py
new file mode 100644
index 0000000000000000000000000000000000000000..576975125b4a3d7433f922896743553beb8e6197
--- /dev/null
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/plotSpectrum.py
@@ -0,0 +1,339 @@
+import numpy as np
+import h5py
+import unyt
+import matplotlib.pyplot as plt
+import argparse
+
+from swiftsimio import load
+from swiftsimio.visualisation.volume_render import render_gas
+
+
+# Parse command line arguments
+argparser = argparse.ArgumentParser()
+argparser.add_argument("input")
+argparser.add_argument("output")
+#argparser.add_argument("resolution")
+args = argparser.parse_args()
+
+# Load snapshot
+filename = args.input
+data = load(filename)
+
+time = np.round(data.metadata.time,2)
+print("Plotting at ", time)
+
+Lbox = data.metadata.boxsize
+
+# Retrieve some information about the simulation run
+artDiffusion = data.metadata.hydro_scheme["Artificial Diffusion Constant"]
+dedHyp = data.metadata.hydro_scheme["Dedner Hyperbolic Constant"]
+dedHypDivv = data.metadata.hydro_scheme["Dedner Hyperbolic div(v) Constant"]
+dedPar = data.metadata.hydro_scheme["Dedner Parabolic Constant"]
+eta = data.metadata.hydro_scheme["Resistive Eta"]
+git = data.metadata.code["Git Revision"]
+gitBranch = data.metadata.code["Git Branch"]
+hydroScheme = data.metadata.hydro_scheme["Scheme"]
+kernel = data.metadata.hydro_scheme["Kernel function"]
+neighbours = data.metadata.hydro_scheme["Kernel target N_ngb"]
+
+# Retrieve particle attributes of interest
+v = data.gas.velocities
+rho = data.gas.densities
+B = data.gas.magnetic_flux_densities
+h = data.gas.smoothing_lengths
+normB = np.sqrt(B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2)
+
+divB = data.gas.magnetic_divergences
+minh = np.min(h.value)
+
+# Renormalize quantities
+mu0 = 1.25663706127e-6 * unyt.kg * unyt.m / (unyt.s ** 2 * unyt.statA ** 2)
+# note that older swiftsimio version is using statA even if swift has MF as A
+#mu0 = 1.25663706127e-1 * unyt.g * unyt.cm / (unyt.s ** 2 * unyt.statA ** 2)
+v = v * (rho[:, None]/2)**0.5
+B = B / np.sqrt(2*mu0)
+
+#Npside = int(len(h)**(1/3))+1
+
+# Generate mass weighted maps of quantities of interest
+data.gas.mass_weighted_vx = data.gas.masses * v[:,0]
+data.gas.mass_weighted_vy = data.gas.masses * v[:,1]
+data.gas.mass_weighted_vz = data.gas.masses * v[:,2]
+
+data.gas.mass_weighted_Bx = data.gas.masses * B[:,0]
+data.gas.mass_weighted_By = data.gas.masses * B[:,1]
+data.gas.mass_weighted_Bz = data.gas.masses * B[:,2]
+
+data.gas.mass_weighted_error = data.gas.masses * np.maximum(h * abs(divB) / (normB + 0.01 * np.max(normB)), 1e-6)
+
+#res = 128 #Npside #int(args.resolution)
+
+#Lslice = 0.5*Lbox
+Lslice_kPc = 20 # 2*21.5
+Lslice = Lslice_kPc * 3.08e18 * 1e3 * unyt.cm
+
+k_slice = 2*np.pi/Lslice
+k_res = np.max(2*np.pi/h)
+
+res = int((2*Lslice/np.min(h)).value) #Npside #int(args.resolution)
+
+# resolution limiter
+resmax= 512
+print('Suggested maximal resoluiton: ',res)
+res = min([res,resmax])
+
+print('Spectral resolution: ',res)
+
+center = Lbox/2
+visualise_region = [
+ center[0] - Lslice,
+ center[0] + Lslice,
+ center[1] - Lslice,
+ center[1] + Lslice,
+ center[2] - Lslice,
+ center[2] + Lslice,
+]
+
+common_arguments = dict(
+ data=data, resolution=res, parallel=True,periodic=True,region=visualise_region,
+)
+
+
+mass_cube = render_gas(**common_arguments, project="masses")
+min_nonzero_m =np.min( mass_cube.flatten()[mass_cube.flatten().value!=0.0])
+mass_cube[mass_cube.value==0]=1e-2*min_nonzero_m
+
+mass_weighted_vx_cube = render_gas(**common_arguments, project="mass_weighted_vx")
+mass_weighted_vy_cube = render_gas(**common_arguments, project="mass_weighted_vy")
+mass_weighted_vz_cube = render_gas(**common_arguments, project="mass_weighted_vz")
+
+mass_weighted_Bx_cube = render_gas(**common_arguments, project="mass_weighted_Bx")
+mass_weighted_By_cube = render_gas(**common_arguments, project="mass_weighted_By")
+mass_weighted_Bz_cube = render_gas(**common_arguments, project="mass_weighted_Bz")
+
+mass_weighted_error_cube = render_gas(**common_arguments, project="mass_weighted_error")
+
+
+vx_cube = mass_weighted_vx_cube/mass_cube
+vy_cube = mass_weighted_vy_cube/mass_cube
+vz_cube = mass_weighted_vz_cube/mass_cube
+
+Bx_cube = mass_weighted_Bx_cube/mass_cube
+By_cube = mass_weighted_By_cube/mass_cube
+Bz_cube = mass_weighted_Bz_cube/mass_cube
+
+error_cube = mass_weighted_error_cube/mass_cube
+
+unit_energy = unyt.g * unyt.cm**2 / unyt.s**2
+unit_energy_density = unit_energy / unyt.cm**3
+unit_sqrt_energy_density = (unit_energy_density)**0.5
+unit_length = 1e3*unyt.pc
+
+vx_cube.convert_to_units(unit_sqrt_energy_density)
+vy_cube.convert_to_units(unit_sqrt_energy_density)
+vz_cube.convert_to_units(unit_sqrt_energy_density)
+
+Bx_cube.convert_to_units(unit_sqrt_energy_density)
+By_cube.convert_to_units(unit_sqrt_energy_density)
+Bz_cube.convert_to_units(unit_sqrt_energy_density)
+
+k_slice = 2*np.pi/Lslice
+k_res = np.max(2*np.pi/h)
+
+k_slice.convert_to_units(1/unit_length)
+k_res.convert_to_units(1/unit_length)
+
+def compute_power_spectrum_scal(Q, dx, nbins):
+ # Grid size
+ Nx, Ny, Nz = Q.shape
+
+ # Compute Fourier transforms
+ Q_k = np.fft.fftn(Q)
+
+ # Compute power spectrum (squared magnitude of Fourier components)
+ Q_power_k = np.abs(Q_k)**2
+
+ # Compute the corresponding wavenumbers
+ kx = np.fft.fftfreq(Nx, d=dx) * 2 * np.pi
+ ky = np.fft.fftfreq(Ny, d=dx) * 2 * np.pi
+ kz = np.fft.fftfreq(Nz, d=dx) * 2 * np.pi
+
+ # Create 3D arrays of wavevectors
+ KX, KY, KZ = np.meshgrid(kx, ky, kz, indexing='ij')
+ k_mag = np.sqrt(KX**2 + KY**2 + KZ**2)
+
+ # Flatten arrays for binning
+ k_mag_flat = k_mag.flatten()
+ Q_power_flat = Q_power_k.flatten()
+
+ # Define k bins (you can tweak bin size)
+ k_bins = np.linspace(0, np.max(2*np.pi/minh), num=nbins)
+
+ # exclude 0th mode:
+ k_bins = k_bins[1:]
+ k_mag_flat = k_mag_flat[1:]
+ Q_power_flat = Q_power_flat[1:]
+
+ k_bin_centers = 0.5 * (k_bins[1:] + k_bins[:-1])
+
+ # Bin the power spectrum
+ power_spectrum, _ = np.histogram(k_mag_flat, bins=k_bins, weights=Q_power_flat)
+ counts, _ = np.histogram(k_mag_flat, bins=k_bins)
+
+ # Avoid division by zero
+ power_spectrum = np.where(counts > 0, power_spectrum / counts, 0)
+
+ return k_bin_centers, power_spectrum
+
+def compute_power_spectrum_vec(Qx, Qy, Qz, dx,nbins):
+ # Ensure all arrays are unyt arrays with same units
+ dx = dx #.to(unyt.pc)
+ volume_element = dx**3 # single cell volume in cm^3
+
+ # Get shape
+ Nx, Ny, Nz = Qx.shape
+
+ # Fourier transform (keep units)
+ Qx_k = np.fft.fftn(Qx.value) * volume_element.value
+ Qy_k = np.fft.fftn(Qy.value) * volume_element.value
+ Qz_k = np.fft.fftn(Qz.value) * volume_element.value
+
+ # Reapply units
+ Qx_k = Qx_k * Qx.units * volume_element.units
+ Qy_k = Qy_k * Qy.units * volume_element.units
+ Qz_k = Qz_k * Qz.units * volume_element.units
+
+ # Power in k-space
+ Pk = (np.abs(Qx_k)**2 + np.abs(Qy_k)**2 + np.abs(Qz_k)**2)#.to(unyt.g/(unyt.s**2)*unyt.cm**5) # or appropriate units
+
+ # Compute wavenumbers
+ kx = np.fft.fftfreq(Nx, d=dx.value) * 2 * np.pi / dx.units #unyt.cm**-1
+ ky = np.fft.fftfreq(Ny, d=dx.value) * 2 * np.pi / dx.units #unyt.cm**-1
+ kz = np.fft.fftfreq(Nz, d=dx.value) * 2 * np.pi / dx.units #unyt.cm**-1
+ KX, KY, KZ = np.meshgrid(kx, ky, kz, indexing='ij')
+ k_mag = np.sqrt(KX**2 + KY**2 + KZ**2)
+
+ # Calculate monopole component
+ Q_dot_k = Qx_k * KX + Qy_k * KY + Qz_k * KZ
+ Qx_k_div = KX * Q_dot_k / k_mag**2
+ Qy_k_div = KY * Q_dot_k / k_mag**2
+ Qz_k_div = KZ * Q_dot_k / k_mag**2
+ Pk_div = (np.abs(Qx_k_div)**2 + np.abs(Qy_k_div)**2 + np.abs(Qz_k_div)**2)#.to(unyt.g/(unyt.s**2)*unyt.cm**5) # or appropriate units
+
+ # Flatten and bin
+ k_flat = k_mag.flatten() #.to('1/cm')
+ Pk_flat = Pk.flatten() #.to(unyt.g/(unyt.s**2)*unyt.cm**5) # this should be correct after unit reduction
+ Pk_div_flat = Pk_div.flatten()
+
+ # Bin in k-space
+ k_min_log = np.log10(np.min(k_flat.value[k_flat.value!=0]))
+ k_max_log = np.log10(np.max(k_flat.value))
+ k_bins = np.logspace(k_min_log,k_max_log, num=nbins) * k_flat.units #*unyt.cm**-1
+
+ # exclude 0th mode:
+ k_bins = k_bins[1:]
+ k_flat = k_flat[1:]
+ Pk_flat = Pk_flat[1:]
+ Pk_div_flat = Pk_div_flat[1:]
+
+ # converting to Energy per unit wavevector
+ Ek_flat = 4*np.pi * k_flat**2 * Pk_flat
+ Ek_div_flat = 4*np.pi * k_flat**2 * Pk_div_flat
+
+ k_bin_centers = 0.5 * (k_bins[1:] + k_bins[:-1])
+
+ power_spectrum = np.zeros(len(k_bin_centers)) * Ek_flat.units
+ power_spectrum_div = np.zeros(len(k_bin_centers)) * Ek_div_flat.units
+ counts = np.zeros(len(k_bin_centers))
+
+ for i in range(len(k_bin_centers)):
+ in_bin = (k_flat >= k_bins[i]) & (k_flat < k_bins[i+1])
+ power_spectrum[i] = Ek_flat[in_bin].sum()
+ power_spectrum_div[i] = Ek_div_flat[in_bin].sum()
+ counts[i] = in_bin.sum()
+
+ # Normalize per mode (optional)
+ power_spectrum = np.where(counts > 0, power_spectrum / counts, 0 * power_spectrum.units)
+ power_spectrum_div = np.where(counts > 0, power_spectrum_div / counts, 0 * power_spectrum_div.units)
+
+ # mask non-zero
+ mask_zeros = power_spectrum==0
+ k_bin_centers = k_bin_centers[~mask_zeros]
+ power_spectrum = power_spectrum[~mask_zeros]
+ power_spectrum_div = power_spectrum_div[~mask_zeros]
+
+ return k_bin_centers, power_spectrum, power_spectrum_div
+
+dx = 2*Lslice.to(unit_length)/(res)
+
+# plot magnetic field spectrum
+ks, Eb, Eb_div = compute_power_spectrum_vec(Bx_cube,By_cube,Bz_cube, dx = dx, nbins=res-1 )
+# plot velocity spectrum
+ks, Ev, _ = compute_power_spectrum_vec(vx_cube,vy_cube,vz_cube, dx = dx, nbins=res-1 )
+# plot divergence error spectrum
+#ks, Perr = compute_power_spectrum_scal(error_cube, dx = dx, nbins=res-1 )
+
+Eb.convert_to_units(unyt.erg*(1e3*unyt.pc))
+Eb_div.convert_to_units(unyt.erg*(1e3*unyt.pc))
+
+Ev.convert_to_units(unyt.erg*(1e3*unyt.pc))
+ks.convert_to_units(1e-3/unyt.pc)
+
+
+fig, ax = plt.subplots(figsize=(10, 6.2))
+
+ax.plot(ks,Ev.value,color='blue',label='$E_v(k)$')
+ax.plot(ks,Eb.value,color='red',linestyle='solid',label='$E_B(k)$')
+ax.plot(ks,Eb_div.value,color='purple',linestyle='solid',label='$E_{B_{mon}}(k)$')
+#ax.plot(ks,Perr,color='purple',label='$P_{R_{0}}(k)$')
+
+# resolution line
+ax.axvline(x=k_res, color='brown',linestyle='solid',label = r'$k_{\mathrm{res}}$')
+
+# self gravity softening line
+ksoft = 2*np.pi / (0.2 * 1e3 * unyt.pc)
+ksoft.convert_to_units(1e-3/unyt.pc)
+ax.axvline(x=ksoft, color='brown',linestyle='dashdot',label = r'$k_{\mathrm{soft}}$')
+
+# plot spectral lines
+ksmock = np.logspace(np.log10(ksoft.value)-1,np.log10(ksoft.value)-0.5,10)*ksoft.units
+p = -5/3
+Eright = 1e57 * unyt.erg*(1e3*unyt.pc)
+kright = ksoft
+kright.convert_to_units(1e-3/unyt.pc)
+Emock = Eright*(ksmock/kright)**(p)
+Emock.convert_to_units(unyt.erg*(1e3*unyt.pc))
+ax.plot(ksmock,Emock,color='black',linestyle='dashed')#,label = '$k^{-5/3}$')
+
+klabel = ksoft/10**0.45
+klabel.convert_to_units(1e-3/unyt.pc)
+Elabel = Eright*(klabel/kright)**(p)
+Elabel.convert_to_units(unyt.erg*(1e3*unyt.pc))
+ax.text(klabel,Elabel,r'$k^{-5/3}$')
+
+# plot spectral lines
+ksmock = np.logspace(np.log10(k_slice.value),np.log10(k_slice.value)+0.5,10)*k_slice.units
+p = 3/2
+Eleft = 1e57 * unyt.erg*(1e3*unyt.pc)
+kleft = k_slice
+kleft.convert_to_units(1e-3/unyt.pc)
+Emock = Eleft*(ksmock/kleft)**(p)
+Emock.convert_to_units(unyt.erg*(1e3*unyt.pc))
+ax.plot(ksmock,Emock,color='black',linestyle='dashed')#,label = '$k^{-5/3}$')
+
+klabel = k_slice/10**0.15
+klabel.convert_to_units(1e-3/unyt.pc)
+Elabel = Eleft*(klabel/kleft)**(p)
+Elabel.convert_to_units(unyt.erg*(1e3*unyt.pc))
+ax.text(klabel,Elabel,r'$k^{3/2}$')
+
+ax.set_yscale('log')
+ax.set_xscale('log')
+ax.set_xlabel(r'$\mathrm{k}$ $[\mathrm{kpc}^{-1}]$', fontsize=30)
+ax.set_ylabel(r'$\mathrm{P}(\mathrm{k})$ $[\mathrm{erg}\cdot\mathrm{kpc}]$', fontsize=30)
+ax.tick_params(axis='both', labelsize=20)
+#ax.grid()
+ax.legend()
+fig.tight_layout()
+plt.savefig(args.output)
+
diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/run.sh b/examples/MHDTests/CoolingHaloWithSpin_MHD/run.sh
index cfbfb42dd8e7aa3e123e13aa97f276c129d16802..c159430636127e86c78c62cf86189b1b76a7cc4e 100755
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/run.sh
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/run.sh
@@ -2,7 +2,13 @@
# Generate the initial conditions if they are not present.
echo "Generating initial conditions for the MHD isothermal potential box example..."
-python3 makeIC_MWsize.py 45000
+python3 makeIC.py 39500
-# Run SWIFT with external potential, SPH and cooling
-../../../swift --external-gravity --self-gravity --hydro --cooling --threads=4 cooling_halo.yml 2>&1 | tee output.log
+# Run SWIFT with external potential, cooling
+#../../../swift --external-gravity --self-gravity --hydro --cooling --stars --star-formation --threads=4 cooling_halo.yml 2>&1 | tee output.log
+
+# Run SWIFT with external potential, cooling, self-gravity and star formation
+#../../../swift --external-gravity --self-gravity --hydro --cooling --stars --star-formation --threads=4 cooling_halo.yml 2>&1 | tee output.log
+
+# Run SWIFT with external potential, cooling, self-gravity, star formation and feedback
+../../../swift --external-gravity --self-gravity --hydro --cooling --stars --star-formation --feedback --sync --threads=4 cooling_halo.yml 2>&1 | tee output.log