From cb952932ae2e30d1d00b1930b940e99e627cf8c3 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Sun, 7 Oct 2018 15:16:53 +0200
Subject: [PATCH] Added a new example: A small cosmological volume with some
cooling on.
---
.../SmallCosmoVolume/small_cosmo_volume.yml | 5 -
examples/SmallCosmoVolume_cooling/README | 14 ++
examples/SmallCosmoVolume_cooling/getIC.sh | 3 +
.../plotTempEvolution.py | 192 ++++++++++++++++++
examples/SmallCosmoVolume_cooling/run.sh | 14 ++
.../small_cosmo_volume.yml | 63 ++++++
src/cooling/const_lambda/cooling_io.h | 5 +-
src/parallel_io.c | 2 +-
src/serial_io.c | 2 +-
src/single_io.c | 2 +-
10 files changed, 292 insertions(+), 10 deletions(-)
create mode 100644 examples/SmallCosmoVolume_cooling/README
create mode 100755 examples/SmallCosmoVolume_cooling/getIC.sh
create mode 100644 examples/SmallCosmoVolume_cooling/plotTempEvolution.py
create mode 100755 examples/SmallCosmoVolume_cooling/run.sh
create mode 100644 examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml
diff --git a/examples/SmallCosmoVolume/small_cosmo_volume.yml b/examples/SmallCosmoVolume/small_cosmo_volume.yml
index 64e20e61d9..353ab469c8 100644
--- a/examples/SmallCosmoVolume/small_cosmo_volume.yml
+++ b/examples/SmallCosmoVolume/small_cosmo_volume.yml
@@ -56,8 +56,3 @@ InitialConditions:
cleanup_velocity_factors: 1
generate_gas_in_ics: 1 # Generate gas particles from the DM-only ICs
cleanup_smoothing_lengths: 1 # Since we generate gas, make use of the (expensive) cleaning-up procedure.
-
-# Constant lambda cooling function
-LambdaCooling:
- lambda_cgs: 1e-22 # Cooling rate (in cgs units)
- cooling_tstep_mult: 1.0 # Dimensionless pre-factor for the time-step condition
diff --git a/examples/SmallCosmoVolume_cooling/README b/examples/SmallCosmoVolume_cooling/README
new file mode 100644
index 0000000000..a0abad5f81
--- /dev/null
+++ b/examples/SmallCosmoVolume_cooling/README
@@ -0,0 +1,14 @@
+Small LCDM cosmological simulation generated by C. Power. Cosmology
+is WMAP9 and the box is 100Mpc/h in size with 64^3 particles.
+We use a softening length of 1/25th of the mean inter-particle separation.
+
+The ICs have been generated to run with Gadget-2 so we need to switch
+on the options to cancel the h-factors and a-factors at reading time.
+We generate gas from the ICs using SWIFT's internal mechanism and set the
+temperature to the expected gas temperature at this redshift.
+
+The 'plotTempEvolution.py' plots the temperature evolution of the gas
+in the simulated volume.
+
+MD5 checksum of the ICs:
+08736c3101fd738e22f5159f78e6022b small_cosmo_volume.hdf5
diff --git a/examples/SmallCosmoVolume_cooling/getIC.sh b/examples/SmallCosmoVolume_cooling/getIC.sh
new file mode 100755
index 0000000000..3b8136cc5a
--- /dev/null
+++ b/examples/SmallCosmoVolume_cooling/getIC.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/small_cosmo_volume.hdf5
+
diff --git a/examples/SmallCosmoVolume_cooling/plotTempEvolution.py b/examples/SmallCosmoVolume_cooling/plotTempEvolution.py
new file mode 100644
index 0000000000..526a417dc6
--- /dev/null
+++ b/examples/SmallCosmoVolume_cooling/plotTempEvolution.py
@@ -0,0 +1,192 @@
+################################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2018 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/>.
+#
+################################################################################
+
+# Computes the temperature evolution of the gas in a cosmological box
+
+# Physical constants needed for internal energy to temperature conversion
+k_in_J_K = 1.38064852e-23
+mH_in_kg = 1.6737236e-27
+
+# Number of snapshots generated
+n_snapshots = 160
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+import os.path
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 9,
+'legend.fontsize': 9,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (3.15,3.15),
+'figure.subplot.left' : 0.14,
+'figure.subplot.right' : 0.99,
+'figure.subplot.bottom' : 0.12,
+'figure.subplot.top' : 0.99,
+'figure.subplot.wspace' : 0.15,
+'figure.subplot.hspace' : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 2.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+# Read the simulation data
+sim = h5py.File("snap_0000.hdf5", "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"][0]
+kernel = sim["/HydroScheme"].attrs["Kernel function"][0]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"][0]
+eta = sim["/HydroScheme"].attrs["Kernel eta"][0]
+alpha = sim["/HydroScheme"].attrs["Alpha viscosity"][0]
+H_mass_fraction = sim["/HydroScheme"].attrs["Hydrogen mass fraction"][0]
+H_transition_temp = sim["/HydroScheme"].attrs["Hydrogen ionization transition temperature"][0]
+T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
+T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
+git = sim["Code"].attrs["Git Revision"]
+cooling_model = sim["/SubgridScheme"].attrs["Cooling Model"]
+
+if cooling_model == "Constant Lambda":
+ Lambda = sim["/SubgridScheme"].attrs["Lambda [cgs]"][0]
+
+# Cosmological parameters
+H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
+gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
+
+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_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
+
+# Primoridal ean molecular weight as a function of temperature
+def mu(T, H_frac=H_mass_fraction, T_trans=H_transition_temp):
+ if T > T_trans:
+ return 4. / (8. - 5. * (1. - H_frac))
+ else:
+ return 4. / (1. + 3. * H_frac)
+
+# Temperature of some primoridal gas with a given internal energy
+def T(u, H_frac=H_mass_fraction, T_trans=H_transition_temp):
+ T_over_mu = (gas_gamma - 1.) * u * mH_in_kg / k_in_J_K
+ ret = np.ones(np.size(u)) * T_trans
+
+ # Enough energy to be ionized?
+ mask_ionized = (T_over_mu > (T_trans+1) / mu(T_trans+1, H_frac, T_trans))
+ if np.sum(mask_ionized) > 0:
+ ret[mask_ionized] = T_over_mu[mask_ionized] * mu(T_trans*10, H_frac, T_trans)
+
+ # Neutral gas?
+ mask_neutral = (T_over_mu < (T_trans-1) / mu((T_trans-1), H_frac, T_trans))
+ if np.sum(mask_neutral) > 0:
+ ret[mask_neutral] = T_over_mu[mask_neutral] * mu(0, H_frac, T_trans)
+
+ return ret
+
+
+z = np.zeros(n_snapshots)
+a = np.zeros(n_snapshots)
+T_mean = np.zeros(n_snapshots)
+T_std = np.zeros(n_snapshots)
+T_log_mean = np.zeros(n_snapshots)
+T_log_std = np.zeros(n_snapshots)
+T_median = np.zeros(n_snapshots)
+T_min = np.zeros(n_snapshots)
+T_max = np.zeros(n_snapshots)
+
+# Loop over all the snapshots
+for i in range(n_snapshots):
+ sim = h5py.File("snap_%04d.hdf5"%i, "r")
+
+ z[i] = sim["/Cosmology"].attrs["Redshift"][0]
+ a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
+
+ u = sim["/PartType0/InternalEnergy"][:]
+
+ # Compute the temperature
+ u *= (unit_length_in_si**2 / unit_time_in_si**2)
+ u /= a[i]**(3 * (gas_gamma - 1.))
+ Temp = T(u)
+
+ # Gather statistics
+ T_median[i] = np.median(Temp)
+ T_mean[i] = Temp.mean()
+ T_std[i] = Temp.std()
+ T_log_mean[i] = np.log10(Temp).mean()
+ T_log_std[i] = np.log10(Temp).std()
+ T_min[i] = Temp.min()
+ T_max[i] = Temp.max()
+
+# CMB evolution
+a_evol = np.logspace(-3, 0, 60)
+T_cmb = (1. / a_evol)**2 * 2.72
+
+# Plot the interesting quantities
+figure()
+subplot(111, xscale="log", yscale="log")
+
+fill_between(a, T_mean-T_std, T_mean+T_std, color='C0', alpha=0.1)
+plot(a, T_max, ls='-.', color='C0', lw=1., label="${\\rm max}~T$")
+plot(a, T_min, ls=':', color='C0', lw=1., label="${\\rm min}~T$")
+plot(a, T_mean, color='C0', label="${\\rm mean}~T$", lw=1.5)
+fill_between(a, 10**(T_log_mean-T_log_std), 10**(T_log_mean+T_log_std), color='C1', alpha=0.1)
+plot(a, 10**T_log_mean, color='C1', label="${\\rm mean}~{\\rm log} T$", lw=1.5)
+plot(a, T_median, color='C2', label="${\\rm median}~T$", lw=1.5)
+
+legend(loc="upper left", frameon=False, handlelength=1.5)
+
+# Cooling model
+if cooling_model == "Constant Lambda":
+ text(1e-2, 6e4, "$\Lambda_{\\rm const} = %.1f\\times10^{%d}~[\\rm{cgs}]$"%(Lambda/10.**(int(log10(Lambda))), log10(Lambda)), fontsize=8)
+else:
+ text(1e-2, 6e4, "No cooling")
+
+# Expected lines
+plot([1e-10, 1e10], [H_transition_temp, H_transition_temp], 'k--', lw=0.5, alpha=0.7)
+text(2.5e-2, H_transition_temp*1.07, "$T_{\\rm HII\\rightarrow HI}$", va="bottom", alpha=0.7, fontsize=8)
+plot([1e-10, 1e10], [T_minimal, T_minimal], 'k--', lw=0.5, alpha=0.7)
+text(1e-2, T_minimal*0.8, "$T_{\\rm min}$", va="top", alpha=0.7, fontsize=8)
+plot(a_evol, T_cmb, 'k--', lw=0.5, alpha=0.7)
+text(a_evol[20], T_cmb[20]*0.55, "$(1+z)^2\\times T_{\\rm CMB,0}$", rotation=-34, alpha=0.7, fontsize=8, va="top", bbox=dict(facecolor='w', edgecolor='none', pad=1.0, alpha=0.9))
+
+
+redshift_ticks = np.array([0., 1., 2., 5., 10., 20., 50., 100.])
+redshift_labels = ["$0$", "$1$", "$2$", "$5$", "$10$", "$20$", "$50$", "$100$"]
+a_ticks = 1. / (redshift_ticks + 1.)
+
+xticks(a_ticks, redshift_labels)
+minorticks_off()
+
+xlabel("${\\rm Redshift}~z$", labelpad=0)
+ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=0)
+xlim(9e-3, 1.1)
+ylim(20, 2.5e7)
+
+savefig("Temperature_evolution.png", dpi=200)
+
diff --git a/examples/SmallCosmoVolume_cooling/run.sh b/examples/SmallCosmoVolume_cooling/run.sh
new file mode 100755
index 0000000000..454ed27950
--- /dev/null
+++ b/examples/SmallCosmoVolume_cooling/run.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e small_cosmo_volume.hdf5 ]
+then
+ echo "Fetching initial conditions for the small cosmological volume example..."
+ ./getIC.sh
+fi
+
+# Run SWIFT
+../swift -c -s -G -C -t 8 small_cosmo_volume.yml 2>&1 | tee output.log
+
+# Plot the temperature evolution
+python plotTempEvolution.py
diff --git a/examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml b/examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml
new file mode 100644
index 0000000000..8f6a2f2ad6
--- /dev/null
+++ b/examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml
@@ -0,0 +1,63 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun
+ UnitLength_in_cgs: 3.08567758e24 # 1 Mpc
+ UnitVelocity_in_cgs: 1e5 # 1 km/s
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+Cosmology: # WMAP9 cosmology
+ Omega_m: 0.276
+ Omega_lambda: 0.724
+ Omega_b: 0.0455
+ h: 0.703
+ a_begin: 0.019607843 # z_ini = 50.
+ a_end: 1.0 # z_end = 0.
+
+# Parameters governing the time integration
+TimeIntegration:
+ dt_min: 1e-6
+ dt_max: 1e-2
+
+# Parameters for the self-gravity scheme
+Gravity:
+ eta: 0.025
+ theta: 0.3
+ comoving_softening: 0.0889 # 1/25th of the mean inter-particle separation: 88.9 kpc
+ max_physical_softening: 0.0889 # 1/25th of the mean inter-particle separation: 88.9 kpc
+ mesh_side_length: 64
+
+# Parameters of the hydro scheme
+SPH:
+ resolution_eta: 1.2348 # "48 Ngb" with the cubic spline kernel
+ CFL_condition: 0.1
+ initial_temperature: 7075. # (1 + z_ini)^2 * 2.72K
+ minimal_temperature: 100.
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: snap
+ delta_time: 1.02
+ scale_factor_first: 0.02
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1.02
+ scale_factor_first: 0.02
+
+Scheduler:
+ max_top_level_cells: 8
+ cell_split_size: 50
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: small_cosmo_volume.hdf5
+ cleanup_h_factors: 1
+ cleanup_velocity_factors: 1
+ generate_gas_in_ics: 1 # Generate gas particles from the DM-only ICs
+ cleanup_smoothing_lengths: 1 # Since we generate gas, make use of the (expensive) cleaning-up procedure.
+
+# Constant lambda cooling function
+LambdaCooling:
+ lambda_cgs: 1e-30 # Cooling rate (in cgs units)
+ cooling_tstep_mult: 1.0 # Dimensionless pre-factor for the time-step condition
diff --git a/src/cooling/const_lambda/cooling_io.h b/src/cooling/const_lambda/cooling_io.h
index 5f517be2a6..e26949e62e 100644
--- a/src/cooling/const_lambda/cooling_io.h
+++ b/src/cooling/const_lambda/cooling_io.h
@@ -35,9 +35,10 @@
* @param h_grpsph The HDF5 group in which to write
*/
__attribute__((always_inline)) INLINE static void cooling_write_flavour(
- hid_t h_grpsph) {
+ hid_t h_grp, const struct cooling_function_data* cooling) {
- io_write_attribute_s(h_grpsph, "Cooling Model", "Constant Lambda");
+ io_write_attribute_s(h_grp, "Cooling Model", "Constant Lambda");
+ io_write_attribute_d(h_grp, "Lambda [cgs]", cooling->lambda_cgs);
}
#endif
diff --git a/src/parallel_io.c b/src/parallel_io.c
index fa1b58d4dc..cb02d1806c 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -1020,7 +1020,7 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
h_grp = H5Gcreate(h_file, "/SubgridScheme", H5P_DEFAULT, H5P_DEFAULT,
H5P_DEFAULT);
if (h_grp < 0) error("Error while creating subgrid group");
- cooling_write_flavour(h_grp);
+ cooling_write_flavour(h_grp, e->cooling_func);
chemistry_write_flavour(h_grp);
H5Gclose(h_grp);
diff --git a/src/serial_io.c b/src/serial_io.c
index 9f874ca8c5..32329a296f 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -870,7 +870,7 @@ void write_output_serial(struct engine* e, const char* baseName,
h_grp = H5Gcreate(h_file, "/SubgridScheme", H5P_DEFAULT, H5P_DEFAULT,
H5P_DEFAULT);
if (h_grp < 0) error("Error while creating subgrid group");
- cooling_write_flavour(h_grp);
+ cooling_write_flavour(h_grp, e->cooling_func);
chemistry_write_flavour(h_grp);
H5Gclose(h_grp);
diff --git a/src/single_io.c b/src/single_io.c
index 33981b6cfa..764d58cf74 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -720,7 +720,7 @@ void write_output_single(struct engine* e, const char* baseName,
h_grp = H5Gcreate(h_file, "/SubgridScheme", H5P_DEFAULT, H5P_DEFAULT,
H5P_DEFAULT);
if (h_grp < 0) error("Error while creating subgrid group");
- cooling_write_flavour(h_grp);
+ cooling_write_flavour(h_grp, e->cooling_func);
chemistry_write_flavour(h_grp);
H5Gclose(h_grp);
--
GitLab