Commit ff8b655c authored by Loic Hausammann's avatar Loic Hausammann Committed by Matthieu Schaller
Browse files

Gear

parent 3a0b20b0
......@@ -1439,13 +1439,14 @@ case "$with_subgrid" in
;;
GEAR)
with_subgrid_cooling=grackle_0
with_subgrid_chemistry=GEAR
with_subgrid_chemistry=GEAR_10
with_subgrid_pressure_floor=GEAR
with_subgrid_stars=GEAR
with_subgrid_star_formation=GEAR
with_subgrid_feedback=none
with_subgrid_feedback=GEAR
with_subgrid_black_holes=none
with_subgrid_task_order=GEAR
# GEAR's order is not used anymore
with_subgrid_task_order=none
enable_fof=no
;;
EAGLE)
......@@ -1791,7 +1792,8 @@ esac
# chemistry function
AC_ARG_WITH([chemistry],
[AS_HELP_STRING([--with-chemistry=<function>],
[chemistry function @<:@none, GEAR, EAGLE default: none@:>@]
[chemistry function @<:@none, GEAR_*, EAGLE default: none@:>@
For GEAR, you need to provide the number of elements (e.g. GEAR_10)]
)],
[with_chemistry="$withval"],
[with_chemistry="none"]
......@@ -1809,8 +1811,10 @@ case "$with_chemistry" in
none)
AC_DEFINE([CHEMISTRY_NONE], [1], [No chemistry function])
;;
GEAR)
GEAR_*)
AC_DEFINE([CHEMISTRY_GEAR], [1], [Chemistry taken from the GEAR model])
number_element=${with_chemistry:5}
AC_DEFINE_UNQUOTED([GEAR_CHEMISTRY_ELEMENT_COUNT], [$number_element], [Number of element to follow])
;;
EAGLE)
AC_DEFINE([CHEMISTRY_EAGLE], [1], [Chemistry taken from the EAGLE model])
......@@ -1913,6 +1917,9 @@ case "$with_feedback" in
EAGLE)
AC_DEFINE([FEEDBACK_EAGLE], [1], [EAGLE stellar feedback and evolution model])
;;
GEAR)
AC_DEFINE([FEEDBACK_GEAR], [1], [GEAR stellar feedback and evolution model])
;;
none)
AC_DEFINE([FEEDBACK_NONE], [1], [No feedback])
;;
......@@ -2139,6 +2146,12 @@ AM_CONDITIONAL([HAVEEAGLECOOLING], [test $with_cooling = "EAGLE"])
# Check if using EAGLE feedback
AM_CONDITIONAL([HAVEEAGLEFEEDBACK], [test $with_feedback = "EAGLE"])
# check if using grackle cooling
AM_CONDITIONAL([HAVEGRACKLECOOLING], [test ${with_cooling:0:7} == "grackle"])
# check if using gear feedback
AM_CONDITIONAL([HAVEGEARFEEDBACK], [test $with_feedback == "GEAR"])
# Handle .in files.
AC_CONFIG_FILES([Makefile src/Makefile examples/Makefile examples/Cooling/CoolingRates/Makefile doc/Makefile doc/Doxyfile tests/Makefile])
AC_CONFIG_FILES([argparse/Makefile tools/Makefile logger/Makefile logger/tests/Makefile])
......
......@@ -7,15 +7,25 @@ InternalUnitSystem:
UnitTemp_in_cgs: 1 # Kelvin
Scheduler:
max_top_level_cells: 8
max_top_level_cells: 16
cell_extra_sparts: 5000 # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation.
cell_max_size: 8000 # (Optional) Maximal number of interactions per task if we force the split (this is the default value).
cell_sub_size_pair_hydro: 2560 # (Optional) Maximal number of hydro-hydro interactions per sub-pair hydro/star task (this is the default value).
cell_sub_size_self_hydro: 3200 # (Optional) Maximal number of hydro-hydro interactions per sub-self hydro/star task (this is the default value).
cell_sub_size_pair_stars: 2560 # (Optional) Maximal number of hydro-star interactions per sub-pair hydro/star task (this is the default value).
cell_sub_size_self_stars: 3200 # (Optional) Maximal number of hydro-star interactions per sub-self hydro/star task (this is the default value).
cell_sub_size_pair_grav: 2560 # (Optional) Maximal number of interactions per sub-pair gravity task (this is the default value).
cell_sub_size_self_grav: 3200 # (Optional) Maximal number of interactions per sub-self gravity task (this is the default value).
cell_split_size: 100 # (Optional) Maximal number of particles per cell (this is the default value).
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 0.5 # The end time of the simulation (in internal units).
dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
dt_max: 0.1 # The maximal time-step size of the simulation (in internal units).
max_dt_RMS_factor: 0.25 # (Optional) Dimensionless factor for the maximal displacement allowed based on the RMS velocities.
# Parameters governing the snapshots
Snapshots:
basename: agora_disk # Common part of the name of output files
......@@ -29,10 +39,12 @@ Statistics:
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
eta: 0.05 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
comoving_softening: 0.08 # Comoving softening length (in internal units).
max_physical_softening: 0.08 # Physical softening length (in internal units).
comoving_DM_softening: 0.08 # Comoving softening length (in internal units).
max_physical_DM_softening: 0.08 # Physical softening length (in internal units).
comoving_baryon_softening: 0.08 # Comoving softening length (in internal units).
max_physical_baryon_softening: 0.08 # Physical softening length (in internal units).
mesh_side_length: 32 # Number of cells along each axis for the periodic gravity mesh.
# Parameters for the hydrodynamics scheme
......@@ -40,6 +52,8 @@ 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.
minimal_temperature: 10. # Kelvin
h_min_ratio: 0.1095 # (Optional) Minimal allowed smoothing length in units of the softening. Defaults to 0 if unspecified.
h_max: 10. # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
# Parameters related to the initial conditions
InitialConditions:
......@@ -55,13 +69,34 @@ LambdaCooling:
# Cooling with Grackle 2.0
GrackleCooling:
CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
WithUVbackground: 1 # Enable or not the UV background
Redshift: 0 # Redshift to use (-1 means time based redshift)
WithMetalCooling: 1 # Enable or not the metal cooling
ProvideVolumetricHeatingRates: 0 # User provide volumetric heating rates
ProvideSpecificHeatingRates: 0 # User provide specific heating rates
SelfShieldingMethod: 0 # Grackle (<= 3) or Gear self shielding method
OutputMode: 1 # Write in output corresponding primordial chemistry mode
MaxSteps: 1000
ConvergenceLimit: 1e-2
cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
with_UV_background: 1 # Enable or not the UV background
redshift: 0 # Redshift to use (-1 means time based redshift)
with_metal_cooling: 1 # Enable or not the metal cooling
provide_volumetric_heating_rates: 0 # User provide volumetric heating rates
provide_specific_heating_rates: 0 # User provide specific heating rates
self_shielding_method: -1 # Grackle (<= 3) or Gear self shielding method
self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
max_steps: 1000
convergence_limit: 1e-2
thermal_time_myr: 5
GEARStarFormation:
star_formation_efficiency: 0.01 # star formation efficiency (c_*)
maximal_temperature: 1e10 # Upper limit to the temperature of a star forming particle
GEARPressureFloor:
jeans_factor: 10
GEARFeedback:
supernovae_energy_erg: 0.1e51
yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5
discrete_yields: 0
GEARChemistry:
initial_metallicity: 1
scale_initial_metallicity: 1
Restarts:
delta_hours: 72 # (Optional) decimal hours between dumps of restart files.
......@@ -14,7 +14,7 @@ copyfile(filename, out)
f = File(out)
for i in range(6):
name = "PartType{}/ElementAbundance".format(i)
name = "PartType{}/ElementAbundances".format(i)
if name in f:
del f[name]
......@@ -23,13 +23,27 @@ for i in range(NPartType):
if name not in f:
continue
grp = f[name + "/SmoothingLength"]
grp = f[name + "/SmoothingLengths"]
grp[:] *= 1.823
# fix issue due to the name of densities
fields = [("Density", "Densities"),
("Entropies", "Entropies"),
("InternalEnergy", "InternalEnergies"),
("SmoothingLength", "SmoothingLengths")]
for field in fields:
if field[1] in f[name] and field[0] not in f[name]:
f[name + "/" + field[0]] = f[name + "/" + field[1]]
cosmo = f["Cosmology"].attrs
head = f["Header"].attrs
head["OmegaLambda"] = cosmo["Omega_lambda"]
head["Omega0"] = cosmo["Omega_b"]
head["HubbleParam"] = cosmo["H0 [internal units]"]
head["Time"] = head["Time"][0]
f.close()
......@@ -6,4 +6,4 @@ if [ "$#" -ne 1 ]; then
exit
fi
wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/$1.hdf5
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraDisk/$1.hdf5
#!/bin/bash
# cleanup work space
rm snapshot_0000.hdf5 snapshot_0500.hdf5
wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/with_cooling/snapshot_0000.hdf5
wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/with_cooling/snapshot_0500.hdf5
# wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/without_cooling/snapshot_0000.hdf5
# wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/without_cooling/snapshot_0500.hdf5
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraDisk/Gear/snapshot_0000.hdf5
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraDisk/Gear/snapshot_0050.hdf5
......@@ -34,7 +34,7 @@ from yt.analysis_modules.halo_analysis.api import *
from yt.data_objects.particle_filters import add_particle_filter
from scipy.stats import kde
from subprocess import call
#mylog.setLevel(1)
# mylog.setLevel(1)
from yt.utilities.math_utils import get_cyl_r, get_cyl_theta, get_cyl_z
......@@ -53,7 +53,7 @@ draw_PDF = 2###2 # 0/1/2/3 = OFF/ON/ON with
draw_pos_vel_PDF = 4##4 # 0/1/2/3/4 = OFF/ON/ON with 1D profile/ON also with 1D dispersion profile/ON also with separate 1D vertical dispersion profiles
draw_star_pos_vel_PDF = 4##4 # 0/1/2/3/4 = OFF/ON/ON with 1D profile/ON also with 1D dispersion profile/ON also with separate 1D vertical dispersion profiles
draw_rad_height_PDF = 2##2 # 0/1/2/3 = OFF/ON/ON with 1D profile/ON with analytic ftn subtracted
draw_metal_PDF = 0##0 # 0/1 = OFF/ON
draw_metal_PDF = 1##1 # 0/1 = OFF/ON
draw_density_DF = 2#2 # 0/1/2 = OFF/ON/ON with difference plot between 1st and 2nd datasets (when 2, dataset_num should be set to 2)
draw_radius_DF = 1#1 # 0/1 = OFF/ON
draw_star_radius_DF = 2#2 # 0/1/2 = OFF/ON/ON with SFR profile and K-S plot (when 2, this automatically turns on draw_radius_DF)
......@@ -110,9 +110,9 @@ marker_names = ['s', 'o', 'p', 'v', '^', '<', '>', 'h', '*']
# [file_location[1]+'GIZMO/snapshot_temp_000', file_location[1]+'GIZMO/snapshot_temp_100']]]
codes = ['SWIFT', 'GEAR']
filenames = [[["./agora_disk_IC.hdf5", "./agora_disk_500Myr.hdf5"], # Sim-noSFF
["./snapshot_0000.hdf5", "./snapshot_0500.hdf5"]],
["./snapshot_0000.hdf5", "./snapshot_0050.hdf5"]],
[["./agora_disk_IC.hdf5", "./agora_disk_500Myr.hdf5"], # Sim-SFF (with star formation and feedback)
["./snapshot_0000.hdf5", "./snapshot_0500.hdf5"]]] # I did not check the order, they can be switched
["./snapshot_0000.hdf5", "./snapshot_0050.hdf5"]]] # I did not check the order, they can be switched
# codes = ["SWIFT"]
# filenames = [[["./agora_disk_0000.hdf5", "./agora_disk_0050.hdf5"]],
......@@ -139,8 +139,8 @@ filenames = [[["./agora_disk_IC.hdf5", "./agora_disk_500Myr.hdf5"], # Sim-noSFF
# filenames = [[[file_location[0]+'GADGET-3/AGORA_ISO_LOW_DRY/snap_iso_dry_000.hdf5', file_location[0]+'GADGET-3/AGORA_ISO_LOW_DRY/snap_iso_dry_010.hdf5']],
# [[file_location[1]+'GADGET-3/AGORA_ISO_LOW_SF_SNII_Thermal_Chevalier_SFT10/snap_iso_sf_000.hdf5', file_location[1]+'GADGET-3/AGORA_ISO_LOW_SF_SNII_Thermal_Chevalier_SFT10/snap_iso_sf_010.hdf5']]]
# codes = ['GEAR']
# filenames = [[['snapshot_0000', 'snapshot_0500']],
# [['snapshot_0000', 'snapshot_0500']]]
# filenames = [[['snapshot_0000.hdf5', 'snapshot_0500.hdf5']],
# [['snapshot_0000.hdf5', 'snapshot_0500']]]
# codes = ['GIZMO']
# filenames = [[[file_location[0]+'GIZMO/snapshot_temp_000', file_location[0]+'GIZMO/snapshot_temp_100']],
# [[file_location[1]+'GIZMO/snapshot_temp_000', file_location[1]+'GIZMO/snapshot_temp_100']]]
......@@ -439,7 +439,7 @@ for time in range(len(times)):
yield line
if draw_PDF >= 1:
fig_PDF += [plt.figure(figsize=(50, 80))]
grid_PDF += [AxesGrid(fig_PDF[time], (0.01,0.01,0.99,0.99), nrows_ncols = (3, int(math.ceil(len(codes)/3.0))), axes_pad = 0.05, add_all = True, share_all = True,
grid_PDF += [AxesGrid(fig_PDF[time], (0.01,0.01,0.99,0.99), nrows_ncols = (2, int(math.ceil(len(codes)/2.0))), axes_pad = 0.05, add_all = True, share_all = True,
label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.05, aspect = False)]
if draw_pos_vel_PDF >= 1:
fig_pos_vel_PDF += [plt.figure(figsize=(50, 80))]
......@@ -469,7 +469,7 @@ for time in range(len(times)):
rad_height_profiles.append([])
if draw_metal_PDF == 1:
fig_metal_PDF += [plt.figure(figsize=(50, 80))]
grid_metal_PDF += [AxesGrid(fig_metal_PDF[time], (0.01,0.01,0.99,0.99), nrows_ncols = (3, int(math.ceil(len(codes)/3.0))), axes_pad = 0.05, add_all = True, share_all = True,
grid_metal_PDF += [AxesGrid(fig_metal_PDF[time], (0.01,0.01,0.99,0.99), nrows_ncols = (2, int(math.ceil(len(codes)/2.0))), axes_pad = 0.05, add_all = True, share_all = True,
label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.05, aspect = False)]
if draw_density_DF >= 1:
density_DF_xs.append([])
......@@ -576,9 +576,9 @@ for time in range(len(times)):
PartType_StarBeforeFiltered_to_use = "PartType4"
if time != 0:
def _FormationTime(field, data):
return pf.arr(data["PartType4", "BirthTime"].d, 'code_time')
return pf.arr(data["PartType4", "BirthTimes"].d, 'code_time')
pf.add_field(("PartType4", FormationTimeType_to_use), function=_FormationTime, particle_type=True, take_log=False, units="code_time")
MassType_to_use = "Masses"
elif codes[code] == "GEAR":
PartType_Gas_to_use = "PartType0"
......@@ -715,6 +715,20 @@ for time in range(len(times)):
def _metallicity_3(field, data):
return data["deposit", PartType_Gas_to_use+"_smoothed_"+MetallicityType_to_use]
pf.add_field(("gas", "metallicity"), function=_metallicity_3, force_override=True, display_name="Metallicity", particle_type=False, take_log=True, units="")
elif codes[code] == 'SWIFT': # "Metals" in SWIFT is 10-species field ([:,9] is the total metal fraction), so requires a change in _vector_fields in frontends/gadget/io.py: added ("Metals", 10)
def _metallicity_2(field, data):
if len(data[PartType_Gas_to_use, "SmoothedElementAbundances"].shape) == 1:
return data[PartType_Gas_to_use, "SmoothedElementAbundances"]
else:
return data[PartType_Gas_to_use, "SmoothedElementAbundances"][:,9].in_units("") # in_units("") turned out to be crucial!; otherwise code_metallicity will be used and it will mess things up
# We are creating ("Gas", "Metallicity") here, different from ("Gas", "metallicity") which is auto-generated by yt but doesn't work properly
pf.add_field((PartType_Gas_to_use, MetallicityType_to_use), function=_metallicity_2, display_name="Metallicity", particle_type=True, take_log=True, units="")
# Also creating smoothed field following an example in yt-project.org/docs/dev/cookbook/calculating_information.html; use hardcoded num_neighbors as in frontends/gadget/fields.py
fn = add_volume_weighted_smoothed_field(PartType_Gas_to_use, "Coordinates", MassType_to_use, "SmoothingLength", "Density", MetallicityType_to_use, pf.field_info, nneighbors=64)
# Alias doesn't work -- e.g. pf.field_info.alias(("gas", "metallicity"), fn[0]) -- probably because pf=GadgetDataset()?, not load()?; so I add and replace existing ("gas", "metallicity")
def _metallicity_3(field, data):
return data["deposit", PartType_Gas_to_use+"_smoothed_"+MetallicityType_to_use]
pf.add_field(("gas", "metallicity"), function=_metallicity_3, force_override=True, display_name="Metallicity", particle_type=False, take_log=True, units="")
if draw_metal_map >= 2:
def _metal_mass(field, data):
return data["gas", "cell_mass"] * data["gas", "metallicity"]
......@@ -745,7 +759,7 @@ for time in range(len(times)):
if draw_metal_map >= 1 or draw_metal_PDF == 1:
def _Metallicity_2(field, data):
return data[(PartType_Gas_to_use, MetallicityType_to_use)]
pf.add_field((PartType_Gas_to_use, "Metallicity_2"), function=_Metallicity_2, take_log=True, particle_type=False, display_name="Metallicity", units="")
pf.add_field((PartType_Gas_to_use, "Metallicity_2"), function=_Metallicity_2, take_log=True, particle_type=True, display_name="Metallicity", units="")
def _Density_2_minus_analytic(field, data):
return data[(PartType_Gas_to_use, "density")] - rho_agora_disk(data[(PartType_Gas_to_use, "radius")], data[(PartType_Gas_to_use, "particle_position_cylindrical_z_abs")])
pf.add_field((PartType_Gas_to_use, "Density_2_minus_analytic"), function=_Density_2_minus_analytic, take_log=False, particle_type=False, display_name="Density Residual", units="g/cm**3")
......@@ -1565,13 +1579,18 @@ for time in range(len(times)):
else:
# Because ParticlePhasePlot doesn't yet work for a log-log PDF for some reason, I will do the following trick.
pf.field_info[(PartType_Gas_to_use, "Metallicity_2")].take_log = False
p3 = PhasePlot(sp, (PartType_Gas_to_use, "density"), (PartType_Gas_to_use, "temperature"), (PartType_Gas_to_use, "Metallicity_2"), weight_field=(PartType_Gas_to_use, MassType_to_use), fontsize=12, x_bins=300, y_bins=300)
p3 = ParticlePhasePlot(sp, (PartType_Gas_to_use, "density"), (PartType_Gas_to_use, "temperature"), (PartType_Gas_to_use, "Metallicity_2"),
weight_field=(PartType_Gas_to_use, MassType_to_use), fontsize=12, x_bins=300, y_bins=300)
p3.set_zlim((PartType_Gas_to_use, "Metallicity_2"), 0.01, 0.04)
p3.set_cmap((PartType_Gas_to_use, "Metallicity_2"), my_cmap2)
plot3 = p3.plots[(PartType_Gas_to_use, "Metallicity_2")]
p3.set_unit("density", "g/cm**3")
p3.set_xlim(1e-29, 1e-21)
p3.set_unit("temperature", "K")
p3.set_ylim(10, 1e7)
p3.set_log("density", True)
p3.set_log("temperature", True)
plot3.figure = fig_metal_PDF[time]
plot3.axes = grid_metal_PDF[time][code].axes
......@@ -2271,7 +2290,7 @@ for time in range(len(times)):
if draw_metal_PDF == 1:
fig_metal_PDF[time].savefig("metal_PDF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
if draw_density_DF >= 1:
plt.clf()
# plt.clf()
# plt.subplot(111, aspect=1)
fig = plt.figure(figsize=(8, 8))
gridspec.GridSpec(4, 1)
......
......@@ -5,6 +5,8 @@
# currently only the low resolution is available
sim=low
rm agora_disk_0*.hdf5
# make run.sh fail if a subcommand fails
set -e
......@@ -19,7 +21,13 @@ fi
if [ ! -e CloudyData_UVB=HM2012.h5 ]
then
echo "Fetching the Cloudy tables required by Grackle..."
../../Cooling/getGrackleCoolingTable.sh
../../Cooling/getGrackleCoolingTable.sh
fi
if [ ! -e chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5 ]
then
echo "Fetching the chemistry tables..."
../getChemistryTable.sh
fi
# copy the initial conditions
......@@ -28,7 +36,7 @@ cp $sim.hdf5 agora_disk.hdf5
python3 changeType.py agora_disk.hdf5
# Run SWIFT
../../swift --cooling --hydro --self-gravity --threads=4 agora_disk.yml 2>&1 | tee output.log
../../swift --sync --limiter --cooling --hydro --self-gravity --star-formation --feedback --stars --threads=8 agora_disk.yml 2>&1 | tee output.log
echo "Changing smoothing length to be Gadget compatible"
......
......@@ -68,3 +68,31 @@ InitialConditions:
cleanup_h_factors: 1 # Remove the h-factors inherited from Gadget
cleanup_velocity_factors: 1 # Remove the sqrt(a) factor in the velocities inherited from Gadget
# Cooling with Grackle 3.0
GrackleCooling:
cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
with_UV_background: 1 # Enable or not the UV background
redshift: -1 # Redshift to use (-1 means time based redshift)
with_metal_cooling: 1 # Enable or not the metal cooling
provide_volumetric_heating_rates: 0 # (optional) User provide volumetric heating rates
provide_specific_heating_rates: 0 # (optional) User provide specific heating rates
self_shielding_method: 0 # (optional) Grackle (<= 3) or Gear self shielding method
max_steps: 10000 # (optional) Max number of step when computing the initial composition
convergence_limit: 1e-2 # (optional) Convergence threshold (relative) for initial composition
thermal_time_Myr: 5
GearChemistry:
initial_metallicity: 0.01295
GEARFeedback:
supernovae_energy_erg: 1e50
yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5
GearPressureFloor:
jeans_factor: 10. # Number of particles required to suppose a resolved clump and avoid the pressure floor.
GEARStarFormation:
star_formation_efficiency: 0.01 # star formation efficiency (c_*)
maximal_temperature: 3e4 # Upper limit to the temperature of a star forming particle
......@@ -7,5 +7,5 @@ then
./getIC.sh
fi
../../swift --feedback --self-gravity --hydro --stars --threads=8 $@ dwarf_galaxy.yml 2>&1 | tee output.log
../../swift --feedback --limiter --sync --self-gravity --hydro --stars --cooling --star-formation --threads=8 $@ dwarf_galaxy.yml 2>&1 | tee output.log
......@@ -7,5 +7,5 @@ then
./getIC.sh
fi
../../swift --feedback --cosmology --self-gravity --hydro --stars --threads=8 zoom_in.yml 2>&1 | tee output.log
../../swift --feedback --cosmology --limiter --sync --self-gravity --hydro --stars --star-formation --threads=8 zoom_in.yml 2>&1 | tee output.log
......@@ -60,3 +60,30 @@ InitialConditions:
cleanup_h_factors: 1 # Remove the h-factors inherited from Gadget
cleanup_velocity_factors: 1 # Remove the sqrt(a) factor in the velocities inherited from Gadget
# Cooling with Grackle 3.0
GrackleCooling:
cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
with_UV_background: 1 # Enable or not the UV background
redshift: -1 # Redshift to use (-1 means time based redshift)
with_metal_cooling: 1 # Enable or not the metal cooling
provide_volumetric_heating_rates: 0 # (optional) User provide volumetric heating rates
provide_specific_heating_rates: 0 # (optional) User provide specific heating rates
self_shielding_method: 0 # (optional) Grackle (<= 3) or Gear self shielding method
max_steps: 10000 # (optional) Max number of step when computing the initial composition
convergence_limit: 1e-2 # (optional) Convergence threshold (relative) for initial composition
thermal_time_Myr: 5
GearChemistry:
initial_metallicity: 0.01295
GEARFeedback:
supernovae_energy_erg: 1e50
yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5
GEARStarFormation:
star_formation_efficiency: 0.01 # star formation efficiency (c_*)
maximal_temperature: 3e4 # Upper limit to the temperature of a star forming particle
GEARPressureFloor:
jeans_factor: 10. # Number of particles required to suppose a resolved clump and avoid the pressure floor.
#!/bin/bash
wget https://obswww.unige.ch/lastro/projects/Clastro/PySWIFTsim/chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1 # Grams
UnitLength_in_cgs: 1 # Centimeters
UnitVelocity_in_cgs: 1 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Values of some physical constants
PhysicalConstants:
G: 0 # (Optional) Overwrite the value of Newton's constant used internally by the code.
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 5e-2 # The end time of the simulation (in internal units).
dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: SN_feedback # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 1e-2 # Time difference between consecutive outputs (in internal units)
compression: 1
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-3 # Time between statistics output
# 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.
# Parameters related to the initial conditions
InitialConditions:
file_name: ./SN_feedback.hdf5
smoothing_length_scaling: 1.
periodic: 1 # Are we running with periodic ICs?
# Parameters for the stellar models
Stars:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_64.hdf5
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 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 h5py
from numpy import *
# Generates a swift IC file for the Sedov blast test in a periodic cubic box
# Parameters
gamma = 5./3. # Gas adiabatic index
rho0 = 1. # Background density
P0 = 1.e-6 # Background pressure
E0= 1. # Energy of the explosion
N_inject = 15 # Number of particles in which to inject energy
fileName = "SN_feedback.hdf5"
#---------------------------------------------------
glass = h5py.File("glassCube_64.hdf5", "r")
# Read particle positions and h from the glass
pos = glass["/PartType0/Coordinates"][:,:]
eps = 1e-6
pos = (pos - pos.min()) / (pos.max() - pos.min() + eps)
h = glass["/PartType0/SmoothingLength"][:] * 0.3 * 3.3
numPart = size(h)
vol = 1.
Boxsize = 1.
# Generate extra arrays
v = zeros((numPart, 3))
ids = linspace(1, numPart, numPart)
m = zeros(numPart)
u = zeros(numPart)
r = zeros(numPart)
r = sqrt((pos[:,0] - 0.5)**2 + (pos[:,1] - 0.5)**2 + (pos[:,2] - 0.5)**2)
m[:] = rho0 * vol / numPart
u[:] = P0 / (rho0 * (gamma - 1))
#--------------------------------------------------
star_pos = zeros((1, 3))
star_pos[:,:] = 0.5 * Boxsize
star_v = zeros((1, 3))
star_v[:,:] = 0.
# increase mass to keep it at center