Skip to content
Snippets Groups Projects
Commit 1317af25 authored by nickishch's avatar nickishch
Browse files

Merge branch... 

Merge branch 'nickishch/karapiperis_antisymmetric_and_faster_Dedner_HaloUpdate' into nickishch/karapiperis_antisymmetric_and_faster_Dedner_OWAR
parents 8694a7f2 ace2185e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 658 additions and 273 deletions
examples/MHDTests/CoolingHaloWithSpin_MHD/1Dprofiles.py
+ 117
− 166
View file @ 1317af25
......@@ -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
......
examples/MHDTests/CoolingHaloWithSpin_MHD/README
+ 10
− 3
View file @ 1317af25
......@@ -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
---------------
......
examples/MHDTests/CoolingHaloWithSpin_MHD/cooling_halo.yml
+ 123
− 70
View file @ 1317af25
# 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.
examples/MHDTests/CoolingHaloWithSpin_MHD/getEaglePhotometryTable.sh 0 → 100755
+ 3
− 0
View file @ 1317af25
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/YieldTables/EAGLE/photometry.tar.gz
tar -xf photometry.tar.gz
examples/MHDTests/CoolingHaloWithSpin_MHD/getYieldTable.sh 0 → 100755
+ 3
− 0
View file @ 1317af25
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/YieldTables/EAGLE/yieldtables.tar.gz
tar -xf yieldtables.tar.gz
examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC.py
+ 49
− 26
View file @ 1317af25
......@@ -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
......
examples/MHDTests/CoolingHaloWithSpin_MHD/plotSolutionReference.py
+ 5
− 5
View file @ 1317af25
......@@ -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),
)
......
examples/MHDTests/CoolingHaloWithSpin_MHD/plotSpectrum.py 0 → 100644
+ 339
− 0
View file @ 1317af25
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)
examples/MHDTests/CoolingHaloWithSpin_MHD/run.sh
+ 9
− 3
View file @ 1317af25
......@@ -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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment