Commit d1824f60 authored by Alexei Borissov's avatar Alexei Borissov

add isolated galaxy feedback example

parent 15790593
Isolated Galaxy generated by the MakeNewDisk code from Springel, Di Matteo &
Hernquist (2005). The done analysis in this example is similar to the work done
by Schaye and Dalla Vecchia (2008) (After this SD08). The default example runs
the simulation for a galaxy with similar mass of their fiducial model and should
produce plots similar to their middle pannel Figure 4. The code needs to be
configured to run with the Hernquist external potential as well as the cooling &
star-formation model of interest. Using the EAGLE model allows to reproduce the
results of SD08.
The code can also be run for other situations to check to verify the law using
different parameters, changes that were done in SD08 are given by:
- gas fraction of 10% instead of 30%, change the IC to f10.hdf5, see getIC.sh,
should reproduce something similar to Figure 4 left hand pannel. Requires
change of fg=.1
- gas fraction of 90% instead of 30%, change the IC to f90.hdf5, see getIC.sh,
should reproduce something similar to Figure 4 right hand pannel. Requires
change of fg=.9
- Changing the effective equation of state to adiabatic, Jeans_gamma_effective
= 1.666667. Should result in something similar to Figure 5 left hand pannel
of SD08.
- Changing the effective equation of state to isothermal, Jeans_gamma_effective
= 1.0000. Should result in something similar to Figure 5 middle hand pannel
of SD08.
- Changing the slope of the Kennicutt-Schmidt law to 1.7, SchmidtLawExponent =
1.7, this should result in a plot similar to Figure 6 of SD08.
- Increasing the density threshold by a factor of 10. thresh_norm_HpCM3 = 1.0,
should reproduce plot similar to Figure 7.
- Decreasing the density threshold by a factor of 10. thresh_norm_HpCM3 = 0.01,
should reproduce plot similar to Figure 7.
- Running with a lower resultion of a factor 8, change the IC to lowres8.hdf5,
see getIC.sh.
- Running with a lower resultion of a factor 64, change the IC to lowres64.hdf5,
see getIC.sh.
- Running with a lower resultion of a factor 512, change the IC to lowres512.hdf5,
see getIC.sh.
Other options to verify the correctness of the code is by chaning the following
parameters:
- Changing the normalization to A/2 or 2A.
- Running the code with zero metallicity.
- Running the code with a factor 6 higher resolution idealized disks, use
highres6.hdf5, see getIC.sh.
- Running with different SPH schemes like Anarchy-PU.
"""
Makes a movie using sphviewer and ffmpeg.
Edit low_frac and up_frac to focus on a certain view of the box.
The colour map can also be changed via colour_map.
Usage: python3 makeMovie.py CoolingHalo_
Written by James Willis (james.s.willis@durham.ac.uk)
"""
import glob
import h5py as h5
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
def getSFH(filename):
# Read the data
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates = f["/PartType4/Coordinates"][:, :]
mass = f["/PartType4/Masses"][:]
# flag = f["/PartType4/NewStarFlag"][:]
birth_time = f["/PartType4/Birth_time"][:]
absmaxz = 2 # kpc
absmaxxy = 10 # kpc
part_mask = (
((coordinates[:, 0] - box_size / 2.0) > -absmaxxy)
& ((coordinates[:, 0] - box_size / 2.0) < absmaxxy)
& ((coordinates[:, 1] - box_size / 2.0) > -absmaxxy)
& ((coordinates[:, 1] - box_size / 2.0) < absmaxxy)
& ((coordinates[:, 2] - box_size / 2.0) > -absmaxz)
& ((coordinates[:, 2] - box_size / 2.0) < absmaxz)
) # & (flag==1)
birth_time = birth_time[part_mask]
mass = mass[part_mask]
histogram = np.histogram(birth_time, bins=200, weights=mass * 2e4, range=[0, 0.1])
values = histogram[0]
xvalues = (histogram[1][:-1] + histogram[1][1:]) / 2.0
return xvalues, values
def getsfrsnapwide():
time = np.arange(1, 101, 1)
SFR_sparticles = np.zeros(100)
SFR_gparticles = np.zeros(100)
new_sparticles = np.zeros(100)
previous_mass = 0
previous_numb = 0
for i in tqdm(range(1, 100)):
# Read the data
filename = "output_%04d.hdf5" % i
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates = f["/PartType0/Coordinates"][:, :]
SFR = f["/PartType0/SFR"][:]
coordinates_star = f["/PartType4/Coordinates"][:, :]
masses_star = f["/PartType4/Masses"][:]
absmaxz = 2 # kpc
absmaxxy = 10 # kpc
part_mask = (
((coordinates[:, 0] - box_size / 2.0) > -absmaxxy)
& ((coordinates[:, 0] - box_size / 2.0) < absmaxxy)
& ((coordinates[:, 1] - box_size / 2.0) > -absmaxxy)
& ((coordinates[:, 1] - box_size / 2.0) < absmaxxy)
& ((coordinates[:, 2] - box_size / 2.0) > -absmaxz)
& ((coordinates[:, 2] - box_size / 2.0) < absmaxz)
& (SFR > 0)
)
SFR = SFR[part_mask]
total_SFR = np.sum(SFR)
SFR_gparticles[i] = total_SFR * 10
return time[:-1], SFR_gparticles[1:]
if __name__ == "__main__":
time, SFR1 = getsfrsnapwide() # , SFR2, SFR_error = getsfrsnapwide()
time2, SFR3 = getSFH("output_%04d.hdf5" % 100)
plt.plot(time2[1:] * 1e3, SFR3[1:], label="Using birth_time of star particles")
plt.plot(time, SFR1, label="Using SFR of gas particles", color="g")
plt.xlabel("Time (Myr)")
plt.ylabel("SFH ($\\rm M_\odot \\rm yr^{-1}$)")
plt.ylim(0, 20)
plt.legend()
plt.savefig("SFH.png")
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/EAGLE/coolingtables.tar.gz
tar -xf coolingtables.tar.gz
#!/bin/bash
#wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/fid.hdf5
# wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/f10.hdf5
# wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/f90.hdf5
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/lowres8.hdf5
# wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/lowres64.hdf5
#wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/lowres512.hdf5
# wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/highres6.hdf5
# Define the system of units to use internally.
InternalUnitSystem:
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: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters for the self-gravity scheme
Gravity:
mesh_side_length: 32 # Number of cells along each axis for the periodic gravity mesh.
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion).
comoving_softening: 0.01 # Comoving softening length (in internal units).
max_physical_softening: 0.01 # Physical softening length (in internal units).
# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 0.1 # The end time of the simulation (in internal units).
dt_min: 1e-9 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: output # Common part of the name of output files
time_first: 0. # (Optional) Time of the first output if non-cosmological time-integration (in internal units)
delta_time: 0.001 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-2 # Time between statistics output
time_first: 0. # (Optional) Time of the first stats output if non-cosmological time-integration (in internal units)
# Parameters related to the initial conditions
InitialConditions:
file_name: lowres8.hdf5 # The file to read
periodic: 0 # Are we running with periodic ICs?
stars_smoothing_length: 0.5
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
h_min_ratio: 0.1 # Minimal smoothing in units of softening.
h_max: 10.
# Standard EAGLE cooling options
EAGLECooling:
dir_name: ./coolingtables/ # Location of the Wiersma+08 cooling tables
H_reion_z: 11.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
# Primordial 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
# Hernquist potential parameters
HernquistPotential:
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)
idealizeddisk: 1 # Run with an idealized galaxy disk
M200: 137.0 # M200 of the galaxy disk
h: 0.704 # reduced Hubble constant (value does not specify the used units!)
concentration: 9.0 # concentration of the Halo
diskfraction: 0.040 # Disk mass fraction
bulgefraction: 0.014 # Bulge mass fraction
timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
epsilon: 0.01 # Softening size (internal units)
# EAGLE star formation parameters
EAGLEStarFormation:
EOS_density_norm_H_p_cm3: 0.1 # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
EOS_temperature_norm_K: 8000 # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
EOS_gamma_effective: 1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
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.
KS_min_over_density: 57.7 # 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.
KS_temperature_margin_dex: 0.5 # Logarithm base 10 of the maximal temperature difference above the EOS allowed to form stars.
threshold_norm_H_p_cm3: 0.1 # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
threshold_Z0: 0.002 # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
threshold_slope: -0.64 # Slope of the metal-dependant star formation threshold
threshold_max_density_H_p_cm3: 10.0 # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
# Parameters for the EAGLE "equation of state"
EAGLEEntropyFloor:
Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
EagleStellarEvolution:
filename: /cosma5/data/Eagle/BG_Tables/YieldTables/
imf_model: Chabrier
EAGLEFeedback:
lifetime_flag: 2
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
# Matthieu Schaller (matthieu.schaller@durham.ac.uk)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
##############################################################################
import matplotlib
matplotlib.use("Agg")
from pylab import *
from scipy import stats
import h5py
import numpy as np
import glob
import os.path
def find_indices(a,b):
result = np.zeros(len(b))
for i in range(len(b)):
result[i] = ((np.where(a == b[i]))[0])[0]
return result
# Plot parameters
params = {'axes.labelsize': 10,
'axes.titlesize': 10,
'font.size': 12,
'legend.fontsize': 12,
'xtick.labelsize': 10,
'ytick.labelsize': 10,
'text.usetex': True,
'figure.figsize' : (9.90,6.45),
'figure.subplot.left' : 0.1,
'figure.subplot.right' : 0.99,
'figure.subplot.bottom' : 0.1,
'figure.subplot.top' : 0.95,
'figure.subplot.wspace' : 0.2,
'figure.subplot.hspace' : 0.2,
'lines.markersize' : 6,
'lines.linewidth' : 3.,
'text.latex.unicode': True
}
rcParams.update(params)
rc('font',**{'family':'sans-serif','sans-serif':['Times']})
# Number of snapshots and elements
newest_snap_name = max(glob.glob('output_*.hdf5'), key=os.path.getctime)
n_snapshots = int(newest_snap_name.replace('output_','').replace('.hdf5','')) + 1
n_elements = 9
# Read the simulation data
sim = h5py.File("output_0000.hdf5", "r")
boxSize = sim["/Header"].attrs["BoxSize"][0]
time = sim["/Header"].attrs["Time"][0]
scheme = sim["/HydroScheme"].attrs["Scheme"]
kernel = sim["/HydroScheme"].attrs["Kernel function"]
neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
eta = sim["/HydroScheme"].attrs["Kernel eta"]
git = sim["Code"].attrs["Git Revision"]
stellar_mass = sim["/PartType4/Masses"][0]
# Units
unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
unit_time_in_cgs = sim["/Units"].attrs["Unit time in cgs (U_t)"]
unit_temp_in_cgs = sim["/Units"].attrs["Unit temperature in cgs (U_T)"]
unit_vel_in_cgs = unit_length_in_cgs / unit_time_in_cgs
unit_energy_in_cgs = unit_mass_in_cgs * unit_vel_in_cgs * unit_vel_in_cgs
unit_length_in_si = 0.01 * unit_length_in_cgs
unit_mass_in_si = 0.001 * unit_mass_in_cgs
unit_time_in_si = unit_time_in_cgs
unit_density_in_cgs = unit_mass_in_cgs*unit_length_in_cgs**-3
unit_pressure_in_cgs = unit_mass_in_cgs/unit_length_in_cgs*unit_time_in_cgs**-2
unit_int_energy_in_cgs = unit_energy_in_cgs/unit_mass_in_cgs
unit_entropy_in_cgs = unit_energy_in_cgs/unit_temp_in_cgs
Gyr_in_cgs = 3.155e16
Msun_in_cgs = 1.989e33
box_energy = zeros(n_snapshots)
box_mass = zeros(n_snapshots)
box_star_mass = zeros(n_snapshots)
box_metal_mass = zeros(n_snapshots)
element_mass = zeros((n_snapshots,n_elements))
t = zeros(n_snapshots)
# Read data from snapshots
for i in range(n_snapshots):
print("reading snapshot "+str(i))
# Read the simulation data
sim = h5py.File("output_%04d.hdf5"%i, "r")
t[i] = sim["/Header"].attrs["Time"][0]
#ids = sim["/PartType0/ParticleIDs"][:]
masses = sim["/PartType0/Masses"][:]
box_mass[i] = np.sum(masses)
star_masses = sim["/PartType4/Masses"][:]
box_star_mass[i] = np.sum(star_masses)
metallicities = sim["/PartType0/Metallicity"][:]
box_metal_mass[i] = np.sum(metallicities * masses)
internal_energies = sim["/PartType0/InternalEnergy"][:]
box_energy[i] = np.sum(masses * internal_energies)
# Plot the interesting quantities
figure()
# Box mass --------------------------------
subplot(221)
plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_mass[1:] - box_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', marker = "*", ms=0.5, label='swift')
xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
ylabel("Change in total gas particle mass (Msun)", labelpad=2)
ticklabel_format(style='sci', axis='y', scilimits=(0,0))
# Box metal mass --------------------------------
subplot(222)
plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_metal_mass[1:] - box_metal_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', ms=0.5, label='swift')
xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
ylabel("Change in total metal mass of gas particles (Msun)", labelpad=2)
ticklabel_format(style='sci', axis='y', scilimits=(0,0))
# Box energy --------------------------------
subplot(223)
plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_energy[1:] - box_energy[0])* unit_energy_in_cgs, linewidth=0.5, color='k', ms=0.5, label='swift')
xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
ylabel("Change in total energy of gas particles (erg)", labelpad=2)
ticklabel_format(style='sci', axis='y', scilimits=(0,0))
# Box mass --------------------------------
subplot(224)
plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_mass[1:] - box_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', marker = "*", ms=0.5, label='gas')
plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_star_mass[1:] - box_star_mass[n_snapshots-1])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='r', marker = "*", ms=0.5, label='stars')
plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_star_mass[1:] - box_star_mass[n_snapshots-1] + box_mass[1:] - box_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='g', marker = "*", ms=0.5, label='total')
xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
ylabel("Change in total gas particle mass (Msun)", labelpad=2)
ticklabel_format(style='sci', axis='y', scilimits=(0,0))
legend()
savefig("box_evolution.png", dpi=200)
#!/bin/bash
if [ ! -e lowres8.hdf5 ]
then
echo "Fetching initial conditions for the isolated galaxy example..."
./getIC.sh
fi
if [ ! -e cooling_tables/ ]
then
echo "Fetching EAGLE cooling tables for the isolated galaxy example..."
./getEagleCoolingTable.sh
fi
../../swift --threads=32 --feedback --external-gravity --self-gravity --stars --star-formation --cooling --temperature --hydro isolated_galaxy.yml 2>&1 | tee output.log
# Kennicutt-Schmidt law plot
python3 plotSolution.py
# Plot that the random star formation matches the expected SFH
python3 SFH.py
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment