diff --git a/examples/MHDTests/CoolingHaloWithSpin_MHD/README b/examples/MHDTests/CoolingHaloWithSpin_MHD/README
index 01b8a8e99d29839deb973e3868012ff60ac642b1..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
+
+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 598284a7e495f733d688eb6f7c3bc4cb7c8d9dc5..c875f3677dc51720e053cc543b74e594799f86e7 100644
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/cooling_halo.yml
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/cooling_halo.yml
@@ -9,7 +9,7 @@ InternalUnitSystem:
# 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).
@@ -20,29 +20,28 @@ Statistics:
recording_triggers_part: [-1, -1] # Not recording as we have many snapshots
recording_triggers_bpart: [-1, -1] # Not recording as we have many snapshots
-#Scheduler:
-# tasks_per_cell: 200
-
# Parameters governing the snapshots
Snapshots:
basename: CoolingHalo # Common part of the name of output files
- time_first: 0. # Time of the first output (in internal units)
- delta_time: 0.1 #0.1 # Time difference between consecutive outputs (in internal units)
-
+ 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 # smoothing length accuracy
- h_min_ratio: 0.1 # Minimal smoothing in units of softening.
h_max: 10.
# Parameters for MHD schemes
MHD:
monopole_subtraction: 1.0
artificial_diffusion: 0.0
- hyperbolic_dedner: 1.0 #1.0
+ hyperbolic_dedner: 1.0
hyperbolic_dedner_divv: 0.5
parabolic_dedner: 1.0
resistive_eta: 0.0
@@ -50,11 +49,9 @@ MHD:
# Parameters related to the initial conditions
InitialConditions:
file_name: CoolingHalo.hdf5 # The file to read
- periodic: 0 # was 1
+ periodic: 0
stars_smoothing_length: 0.5
-# Note: softening 2.6 below is for 1.4e7 Msol particle mass to impose density cut 1e2 cm^-3
-
# External potential parameters
NFWPotential:
useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
@@ -63,16 +60,16 @@ NFWPotential:
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.2 # Softening size (internal units)
+ epsilon: 0.7 # Softening size (internal units)
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
- MAC: geometric
- theta_cr: 0.7 # Opening angle (Multipole acceptance criterion).
- use_tree_below_softening: 1
- max_physical_DM_softening: 0.2 # Physical softening length (in internal units).
- max_physical_baryon_softening: 0.2 # Physical softening length (in internal units).
+ 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:
@@ -134,3 +131,47 @@ EAGLEEntropyFloor:
Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
Cool_temperature_norm_K: 10. # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin. (NOTE: This is below the min T and hence this floor does nothing)
Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
+
+# 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/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 af63329d0a35d9d7e1a327c14f52876ac893ec66..a06466e55816062f2c98cc4ca588f99d60adb6ee 100644
--- a/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC.py
+++ b/examples/MHDTests/CoolingHaloWithSpin_MHD/makeIC.py
@@ -51,7 +51,7 @@ f_b = 0.17 # baryon fraction
c_200 = 7.2 # concentration parameter
# Set the magnitude of the uniform seed magnetic field
-B0_Gaussian_Units = 0.0 #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
@@ -118,9 +118,11 @@ print("G=", const_G)
# Parameters
periodic = 1 # 1 For periodic box
-boxSize = 2.0 # 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
@@ -196,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/3 #boxSize * np.sqrt(3)
+coords *= 2*R_max
# calculate max distance from the center (units of r_200)
-R_max = boxSize/6 #boxSize/2 #np.sqrt(3) * (boxSize / 2)
r_s = 1 / c_200
# calculate distances to the center
r_uni = np.linalg.norm(coords, axis=1)
@@ -465,11 +466,6 @@ elif spin_lambda_choice == "Bullock":
jsp_B = (np.sqrt(2) * v_200 )*spin_lambda
v *= jsp_B/jsp
- # jp_sp = Lp_tot_abs / (gas_particle_mass * np.sum(mask_r_200))
- # calc_spin_par_B = jp_sp / (np.sqrt(2) * 1 * v_200)
- # v *= spin_lambda / calc_spin_par_B
- # l *= spin_lambda / calc_spin_par_B
-
# Draw velocity distribution
if plot_v_distribution:
import matplotlib.pyplot as plt
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/run.sh b/examples/MHDTests/CoolingHaloWithSpin_MHD/run.sh
index 7d3d0c5d24ed0a4f0a02b6f0db5dbedcca9ef912..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 --stars --star-formation --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