diff --git a/Makefile.am b/Makefile.am
index 5b23ff459ad3edda975b3280efe2dbed8f566e04..3fb1a51af67e3c014cbeb47e290df9e75a7ff032 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -21,7 +21,6 @@ ACLOCAL_AMFLAGS = -I m4
# Show the way...
SUBDIRS = src examples doc tests tools
if HAVEEAGLECOOLING
-SUBDIRS += examples/CoolingBox
SUBDIRS += examples/CoolingRates
endif
diff --git a/examples/CoolingBox/Makefile.am b/examples/CoolingBox/Makefile.am
deleted file mode 100644
index ca9ddd8670957273b22cd8d4c6d91db782009f0e..0000000000000000000000000000000000000000
--- a/examples/CoolingBox/Makefile.am
+++ /dev/null
@@ -1,32 +0,0 @@
-# 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 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 General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-# Add the source directory and the non-standard paths to the included library headers to CFLAGS
-AM_CFLAGS = -I$(top_srcdir)/src $(HDF5_CPPFLAGS) $(GSL_INCS) $(FFTW_INCS)
-
-AM_LDFLAGS = $(HDF5_LDFLAGS) $(HDF5_LIBS) $(FFTW_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(TBBMALLOC_LIBS) $(GRACKLE_LIBS) $(GSL_LIBS) $(PROFILER_LIBS)
-
-# Extra libraries.
-EXTRA_LIBS = $(HDF5_LIBS) $(FFTW_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(TBBMALLOC_LIBS) $(GRACKLE_LIBS) $(VELOCIRAPTOR_LIBS) $(GSL_LIBS)
-
-# Programs.
-bin_PROGRAMS = testCooling
-
-# Sources
-testCooling_SOURCES = testCooling.c
-testCooling_CFLAGS = $(MYFLAGS) $(AM_CFLAGS)
-testCooling_LDADD = ../../src/.libs/libswiftsim.a $(EXTRA_LIBS)
-
diff --git a/examples/CoolingBox/README b/examples/CoolingBox/README
deleted file mode 100644
index 1e41bad4afaf690f41a8c75f05b0452fea51b022..0000000000000000000000000000000000000000
--- a/examples/CoolingBox/README
+++ /dev/null
@@ -1,5 +0,0 @@
-This is a test of cooling on a box of uniform gas. It can be run with a variety of of cooling models, both with and without cosmology. If running without cosmology specify coolingBox_non-cosmo.yml
-
-In the case of EAGLE cooling a script (test_script.sh) is provided to run with a variety of parameters and the solution verified against a subcycled explicit solution calculated by the python script analytical_test.py. The parameters that are tested for are: redshift, metal abundance (solar or primordial), hydrogen number density and initial energy. These parameters may be changed in test_script.sh by modifying the for loops, or specifying the desired metal abundances. Currently the script is designed to work with cosmology only. Please note that in order for the cooling box solution to be verified against the subcycled solution, the cooling rates for the parameters specified in test_script.sh are calculated by ../CoolingRates/testCooling.
-
-Initial conditions are produced from a glass file, specified in run.sh and makeIC.py.
diff --git a/examples/CoolingBox/analytical_test.py b/examples/CoolingBox/analytical_test.py
deleted file mode 100644
index c886c58d032b662724df4da81f567705033c3f15..0000000000000000000000000000000000000000
--- a/examples/CoolingBox/analytical_test.py
+++ /dev/null
@@ -1,153 +0,0 @@
-import sys
-import matplotlib
-matplotlib.use("Agg")
-from pylab import *
-import h5py
-
-# 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' : (3.15,3.15),
-'figure.subplot.left' : 0.2,
-'figure.subplot.right' : 0.99,
-'figure.subplot.bottom' : 0.2,
-'figure.subplot.top' : 0.9,
-'figure.subplot.wspace' : 0.15,
-'figure.subplot.hspace' : 0.12,
-'lines.markersize' : 6,
-'lines.linewidth' : 3.,
-'text.latex.unicode': True
-}
-rcParams.update(params)
-rc('font',**{'family':'sans-serif','sans-serif':['Times']})
-
-
-import numpy as np
-import h5py as h5
-
-# function for interpolating 1d table of cooling rates
-def interpol_lambda(u_list,cooling_rate,u):
- if u < u_list[0]:
- return[cooling_rate[0]]
- if u > u_list[len(u_list)-1]:
- return[cooling_rate[len(cooling_rate)-1]]
- j = 0
- while u_list[j+1] < u:
- j += 1
- cooling = cooling_rate[j]*(u_list[j+1]-u)/(u_list[j+1]-u_list[j]) + cooling_rate[j+1]*(u-u_list[j])/(u_list[j+1]-u_list[j])
- return cooling
-
-# File containing the total energy
-stats_filename = "./energy.txt"
-
-# First snapshot
-snap_filename = "coolingBox_0000.hdf5"
-
-# Some constants in cgs units
-k_b = 1.38E-16 #boltzmann
-m_p = 1.67e-24 #proton mass
-
-# Initial conditions set in makeIC.py
-T_init = 1.0e7
-
-# Read the initial state of the gas
-f = h5.File(snap_filename,'r')
-#rho = np.mean(f["/PartType0/Density"])
-#pressure = np.mean(f["/PartType0/Pressure"])
-
-# Read the units parameters from the snapshot
-units = f["InternalCodeUnits"]
-unit_mass = units.attrs["Unit mass in cgs (U_M)"]
-unit_length = units.attrs["Unit length in cgs (U_L)"]
-unit_time = units.attrs["Unit time in cgs (U_t)"]
-unit_vel = unit_length/unit_time
-
-
-# Read snapshots
-if len(sys.argv) >= 4:
- nsnap = int(sys.argv[5])
-else:
- print("Need to specify number of snapshots, defaulting to 100.")
- nsnap = 100
-npart = 4096
-u_snapshots_cgs = zeros(nsnap)
-u_part_snapshots_cgs = zeros((nsnap,npart))
-t_snapshots_cgs = zeros(nsnap)
-scale_factor = zeros(nsnap)
-rho = zeros(nsnap)
-for i in range(nsnap):
- snap = h5.File("coolingBox_%0.4d.hdf5"%i,'r')
- u_part_snapshots_cgs[i][:] = snap["/PartType0/InternalEnergy"][:]*(unit_length**2)/(unit_time**2)
- t_snapshots_cgs[i] = snap["/Header"].attrs["Time"] * unit_time
- scale_factor[i] = snap["/Header"].attrs["Scale-factor"]
- rho[i] = np.mean(snap["/PartType0/Density"])*unit_mass/(unit_length**3)/(scale_factor[i]**3)
-
-
-# calculate analytic solution. Since cooling rate is constant,
-# temperature decrease in linear in time.
-# read Lambda and temperatures from table
-internal_energy = []
-cooling_rate = []
-file_in = open('cooling_output.dat', 'r')
-for line in file_in:
- data = line.split()
- internal_energy.append(float(data[0]))
- cooling_rate.append(-float(data[1]))
-
-tfinal = t_snapshots_cgs[nsnap-1]
-tfirst = t_snapshots_cgs[0]
-nt = nsnap
-dt = (tfinal-tfirst)/nt
-
-t_sol = np.zeros(int(nt))
-u_sol = np.zeros(int(nt))
-lambda_sol = np.zeros(int(nt))
-u = np.mean(u_part_snapshots_cgs[0,:])/scale_factor[0]**2
-u_sol[0] = u
-t_sol[0] = tfirst
-# integrate explicit ODE
-for j in range(int(nt-1)):
- t_sol[j+1] = t_sol[j] + dt
- Lambda_net = interpol_lambda(internal_energy,cooling_rate,u_sol[j])
- if int(sys.argv[4]) == 1:
- nH = 0.702*rho[j]/m_p
- ratefact = nH*0.702/m_p
- else:
- nH = 0.752*rho[j]/m_p
- ratefact = nH*0.752/m_p
- print(rho[j],nH,ratefact,u,Lambda_net)
- u_next = u - Lambda_net*ratefact*dt
- print(u_next, u, ratefact, Lambda_net, dt, nH)
- u_sol[j+1] = u_next
- u = u_next
-u_sol = u_sol*scale_factor[0]**2
-
-# swift solution
-mean_u = np.zeros(nsnap)
-for j in range(nsnap):
- mean_u[j] = np.mean(u_part_snapshots_cgs[j,:])
-
-# plot and save results
-log_nh = float(sys.argv[2])
-redshift = float(sys.argv[1])
-p = plt.figure(figsize = (6,6))
-p1, = plt.loglog(t_sol,u_sol,linewidth = 0.5,color = 'k', marker = '*', markersize = 1.5,label = 'explicit ODE')
-p2, = plt.loglog(t_snapshots_cgs,mean_u,linewidth = 0.5,color = 'r', marker = '*', markersize = 1.5,label = 'swift mean')
-l = legend(handles = [p1,p2])
-xlabel("${\\rm{Time~[s]}}$", labelpad=0, fontsize = 14)
-ylabel("Internal energy ${\\rm{[erg \cdot g^{-1}]}}$", fontsize = 14)
-if int(sys.argv[4]) == 1:
- title('$n_h = 10^{' + "{}".format(log_nh) + '} \\rm{cm}^{-3}, z = ' + "{}".format(redshift) + '$, solar metallicity,\n relative error: ' + "{0:.4f}".format( (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1])), fontsize = 16)
- name = "z_"+str(sys.argv[1])+"_nh_"+str(sys.argv[2])+"_pressure_"+str(sys.argv[3])+"_solar.png"
-elif int(sys.argv[4]) == 0:
- title('$n_h = 10^{' + "{}".format(log_nh) + '} \\rm{cm}^{-3}, z = ' + "{}".format(redshift) + '$, zero metallicity,\n relative error: ' + "{0:.4f}".format( (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1])), fontsize = 16)
- name = "z_"+str(sys.argv[1])+"_nh_"+str(sys.argv[2])+"_pressure_"+str(sys.argv[3])+"_zero_metal.png"
-
-savefig(name, dpi=200)
-
-print('Final internal energy ode, swift, error ' + str(u_sol[int(nt)-1]) + ' ' + str(mean_u[nsnap-1]) + ' ' + str( (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1])))
diff --git a/examples/CoolingBox/coolingBox.yml b/examples/CoolingBox/coolingBox.yml
index 781919a498e1cac7626c42d33352e0c507133a0a..df2c29c0b612eff377423b7bb76e2c8e1e530df1 100644
--- a/examples/CoolingBox/coolingBox.yml
+++ b/examples/CoolingBox/coolingBox.yml
@@ -1,51 +1,28 @@
# Define the system of units to use internally.
InternalUnitSystem:
- UnitMass_in_cgs: 1.989e43 # 10^10 M_sun in grams
- UnitLength_in_cgs: 3.085678e24 # Mpc in centimeters
- UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
-
-# Cosmological parameters
-Cosmology:
- h: 0.6777 # Reduced Hubble constant
- a_begin: 0.04 # Initial scale-factor of the simulation
- a_end: 1.0 # Final scale factor of the simulation
- Omega_m: 0.307 # Matter density parameter
- Omega_lambda: 0.693 # Dark-energy density parameter
- Omega_b: 0.0455 # Baryon density parameter
+ UnitMass_in_cgs: 2.0e33 # Solar masses
+ UnitLength_in_cgs: 3.0857e21 # Kiloparsecs
+ UnitVelocity_in_cgs: 1.0e5 # Kilometers per second
+ UnitCurrent_in_cgs: 1 # 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: 1e-2 # 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-6 # The maximal time-step size of the simulation (in internal units).
-
-Scheduler:
- max_top_level_cells: 15
-
+ time_end: 0.25 # The end time of the simulation (in internal units).
+ dt_min: 1e-5 # 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: coolingBox # Common part of the name of output files
- scale_factor_first: 0.142857142857 # Scale-factor of the first snaphot (cosmological run)
- time_first: 0.01 # Time of the first output (non-cosmological run) (in internal units)
- delta_time: 1.00002 # Time difference between consecutive outputs (in internal units)
- compression: 1
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 1e-2 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
- scale_factor_first: 0.142857142857 # Scale-factor of the first stat dump (cosmological run)
- time_first: 0.01 # Time of the first stat dump (non-cosmological run) (in internal units)
- delta_time: 1.00002 # Time between statistics output
+ delta_time: 1e-3 # Time between statistics output
-# Parameters for the self-gravity scheme
-Gravity:
- eta: 0.025 # Constant dimensionless multiplier for time integration.
- theta: 0.85 # Opening angle (Multipole acceptance criterion)
- comoving_softening: 0.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
-
# 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).
@@ -62,12 +39,6 @@ LambdaCooling:
lambda_nH2_cgs: 1e-22 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3])
cooling_tstep_mult: 1.0 # Dimensionless pre-factor for the time-step condition
-# Dimensionless constant cooling (AB 13/02/18)
-ConstCooling:
- cooling_rate: 10000.0
- min_energy: 0.0
- cooling_tstep_mult: 1.0
-
# Cooling with Grackle 2.0
GrackleCooling:
CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
@@ -80,14 +51,13 @@ GrackleCooling:
OutputMode: 1 # Write in output corresponding primordial chemistry mode
MaxSteps: 1000
ConvergenceLimit: 1e-2
-
+
EagleCooling:
filename: /cosma5/data/Eagle/BG_Tables/CoolingTables/
reionisation_redshift: 8.989
he_reion_z_center: 3.5
he_reion_z_sigma: 0.5
he_reion_ev_pH: 2.0
- newton_integration: 1
EAGLEChemistry:
InitMetallicity: 0.014
@@ -105,4 +75,3 @@ EAGLEChemistry:
GearChemistry:
InitialMetallicity: 0.01295
-
diff --git a/examples/CoolingBox/coolingBox_non-cosmo.yml b/examples/CoolingBox/coolingBox_non-cosmo.yml
deleted file mode 100644
index 2c10909022f0291bf60b4424a53ed27efa47d98a..0000000000000000000000000000000000000000
--- a/examples/CoolingBox/coolingBox_non-cosmo.yml
+++ /dev/null
@@ -1,81 +0,0 @@
-# Define the system of units to use internally.
-InternalUnitSystem:
- UnitMass_in_cgs: 2.0e33 # Solar masses
- UnitLength_in_cgs: 3.0857e21 # Kiloparsecs
- UnitVelocity_in_cgs: 1.0e5 # Kilometers per second
- UnitCurrent_in_cgs: 1 # 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: 0.01 # 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-5 # The maximal time-step size of the simulation (in internal units).
-
-# Parameters governing the snapshots
-Snapshots:
- basename: coolingBox # Common part of the name of output files
- time_first: 0. # Time of the first output (in internal units)
- delta_time: 2e-5 # Time difference between consecutive outputs (in internal units)
-
-# 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: ./coolingBox.hdf5 # The file to read
- periodic: 1
-
-# Dimensionless pre-factor for the time-step condition
-LambdaCooling:
- lambda_cgs: 1.0e-22 # Cooling rate (in cgs units)
- minimum_temperature: 1.0e4 # Minimal temperature (Kelvin)
- mean_molecular_weight: 0.59 # Mean molecular weight
- hydrogen_mass_abundance: 0.75 # Hydrogen mass abundance (dimensionless)
- cooling_tstep_mult: 1.0 # Dimensionless pre-factor for the time-step condition
-
-# Cooling with Grackle 2.0
-GrackleCooling:
- CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
- WithUVbackground: 0 # 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
-
-EagleCooling:
- filename: /cosma5/data/Eagle/BG_Tables/CoolingTables/
- reionisation_redshift: 8.898
- he_reion: 1
- he_reion_z_center: 3.5
- he_reion_z_sigma: 0.5
- he_reion_ev_pH: 2.0
-
-EAGLEChemistry:
- InitMetallicity: 0.014
- InitAbundance_Hydrogen: 0.70649785
- InitAbundance_Helium: 0.28055534
- InitAbundance_Carbon: 2.0665436e-3
- InitAbundance_Nitrogen: 8.3562563e-4
- InitAbundance_Oxygen: 5.4926244e-3
- InitAbundance_Neon: 1.4144605e-3
- InitAbundance_Magnesium: 5.907064e-4
- InitAbundance_Silicon: 6.825874e-4
- InitAbundance_Iron: 1.1032152e-3
- CalciumOverSilicon: 0.0941736
- SulphurOverSilicon: 0.6054160
-
-GearChemistry:
- InitialMetallicity: 0.01295
-
diff --git a/examples/CoolingBox/coolingBox_primordial.yml b/examples/CoolingBox/coolingBox_primordial.yml
deleted file mode 100644
index e352dbdb9f1469dfd939c2ff5fb1431f68b9cc5f..0000000000000000000000000000000000000000
--- a/examples/CoolingBox/coolingBox_primordial.yml
+++ /dev/null
@@ -1,107 +0,0 @@
-# Define the system of units to use internally.
-InternalUnitSystem:
- UnitMass_in_cgs: 1.989e43 # 10^10 M_sun in grams
- UnitLength_in_cgs: 3.085678e24 # Mpc in centimeters
- UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
-
-# Cosmological parameters
-Cosmology:
- h: 0.6777 # Reduced Hubble constant
- a_begin: 0.142857142857 # Initial scale-factor of the simulation
- a_end: 0.144285714286 # Final scale factor of the simulation
- Omega_m: 0.307 # Matter density parameter
- Omega_lambda: 0.693 # Dark-energy density parameter
- Omega_b: 0.0455 # Baryon density parameter
-
-# Parameters governing the time integration
-TimeIntegration:
- time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 1e-2 # 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-9 # The maximal time-step size of the simulation (in internal units).
-
-Scheduler:
- max_top_level_cells: 15
-
-# Parameters governing the snapshots
-Snapshots:
- basename: coolingBox # Common part of the name of output files
- scale_factor_first: 0.142857142857 # Scale-factor of the first snaphot (cosmological run)
- time_first: 0.01 # Time of the first output (non-cosmological run) (in internal units)
- delta_time: 1.00002 # Time difference between consecutive outputs (in internal units)
- compression: 1
-
-# Parameters governing the conserved quantities statistics
-Statistics:
- scale_factor_first: 0.142857142857 # Scale-factor of the first stat dump (cosmological run)
- time_first: 0.01 # Time of the first stat dump (non-cosmological run) (in internal units)
- delta_time: 1.00002 # Time between statistics output
-
-# Parameters for the self-gravity scheme
-Gravity:
- eta: 0.025 # Constant dimensionless multiplier for time integration.
- theta: 0.85 # Opening angle (Multipole acceptance criterion)
- comoving_softening: 0.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
-
-# 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.
- minimal_temperature: 100. # Kelvin
-
-# Parameters related to the initial conditions
-InitialConditions:
- file_name: ./coolingBox.hdf5 # The file to read
- periodic: 1
-
-# Dimensionless pre-factor for the time-step condition
-LambdaCooling:
- lambda_nH2_cgs: 1e-22 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3])
- cooling_tstep_mult: 1.0 # Dimensionless pre-factor for the time-step condition
-
-# Dimensionless constant cooling (AB 13/02/18)
-ConstCooling:
- cooling_rate: 10000.0
- min_energy: 0.0
- cooling_tstep_mult: 1.0
-
-# Cooling with Grackle 2.0
-GrackleCooling:
- CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
- WithUVbackground: 0 # 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
-
-EagleCooling:
- filename: /cosma5/data/Eagle/BG_Tables/CoolingTables/
- reionisation_redshift: 8.989
- he_reion_z_center: 3.5
- he_reion_z_sigma: 0.5
- he_reion_ev_pH: 2.0
-
-EAGLEChemistry:
- InitMetallicity: 0.0
- InitAbundance_Hydrogen: 0.752
- InitAbundance_Helium: 0.248
- InitAbundance_Carbon: 0.000
- InitAbundance_Nitrogen: 0.000
- InitAbundance_Oxygen: 0.000
- InitAbundance_Neon: 0.000
- InitAbundance_Magnesium: 0.000
- InitAbundance_Silicon: 0.000
- InitAbundance_Iron: 0.000
- CalciumOverSilicon: 0.0941736
- SulphurOverSilicon: 0.6054160
-
-GearChemistry:
- InitialMetallicity: 0.01295
-
diff --git a/examples/CoolingBox/coolingBox_solar.yml b/examples/CoolingBox/coolingBox_solar.yml
deleted file mode 100644
index ccdb85eb802d6ee8e1368529839e058a47d165a3..0000000000000000000000000000000000000000
--- a/examples/CoolingBox/coolingBox_solar.yml
+++ /dev/null
@@ -1,107 +0,0 @@
-# Define the system of units to use internally.
-InternalUnitSystem:
- UnitMass_in_cgs: 1.989e43 # 10^10 M_sun in grams
- UnitLength_in_cgs: 3.085678e24 # Mpc in centimeters
- UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
-
-# Cosmological parameters
-Cosmology:
- h: 0.6777 # Reduced Hubble constant
- a_begin: 0.1 # Initial scale-factor of the simulation
- a_end: 1.0 # Final scale factor of the simulation
- Omega_m: 0.307 # Matter density parameter
- Omega_lambda: 0.693 # Dark-energy density parameter
- Omega_b: 0.0455 # Baryon density parameter
-
-# Parameters governing the time integration
-TimeIntegration:
- time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 1e-2 # 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-7 # The maximal time-step size of the simulation (in internal units).
-
-Scheduler:
- max_top_level_cells: 15
-
-# Parameters governing the snapshots
-Snapshots:
- basename: coolingBox # Common part of the name of output files
- scale_factor_first: 0.142857142857 # Scale-factor of the first snaphot (cosmological run)
- time_first: 0.01 # Time of the first output (non-cosmological run) (in internal units)
- delta_time: 1.00002 # Time difference between consecutive outputs (in internal units)
- compression: 1
-
-# Parameters governing the conserved quantities statistics
-Statistics:
- scale_factor_first: 0.142857142857 # Scale-factor of the first stat dump (cosmological run)
- time_first: 0.01 # Time of the first stat dump (non-cosmological run) (in internal units)
- delta_time: 1.00002 # Time between statistics output
-
-# Parameters for the self-gravity scheme
-Gravity:
- eta: 0.025 # Constant dimensionless multiplier for time integration.
- theta: 0.85 # Opening angle (Multipole acceptance criterion)
- comoving_softening: 0.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
-
-# 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.
- minimal_temperature: 100. # Kelvin
-
-# Parameters related to the initial conditions
-InitialConditions:
- file_name: ./coolingBox.hdf5 # The file to read
- periodic: 1
-
-# Dimensionless pre-factor for the time-step condition
-LambdaCooling:
- lambda_nH2_cgs: 1e-22 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3])
- cooling_tstep_mult: 1.0 # Dimensionless pre-factor for the time-step condition
-
-# Dimensionless constant cooling (AB 13/02/18)
-ConstCooling:
- cooling_rate: 10000.0
- min_energy: 0.0
- cooling_tstep_mult: 1.0
-
-# Cooling with Grackle 2.0
-GrackleCooling:
- CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
- WithUVbackground: 0 # 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
-
-EagleCooling:
- filename: /cosma5/data/Eagle/BG_Tables/CoolingTables/
- reionisation_redshift: 8.989
- he_reion_z_center: 3.5
- he_reion_z_sigma: 0.5
- he_reion_ev_pH: 2.0
-
-EAGLEChemistry:
- InitMetallicity: 0.014
- InitAbundance_Hydrogen: 0.70649785
- InitAbundance_Helium: 0.28055534
- InitAbundance_Carbon: 2.0665436e-3
- InitAbundance_Nitrogen: 8.3562563e-4
- InitAbundance_Oxygen: 5.4926244e-3
- InitAbundance_Neon: 1.4144605e-3
- InitAbundance_Magnesium: 5.907064e-4
- InitAbundance_Silicon: 6.825874e-4
- InitAbundance_Iron: 1.1032152e-3
- CalciumOverSilicon: 0.0941736
- SulphurOverSilicon: 0.6054160
-
-GearChemistry:
- InitialMetallicity: 0.01295
-
diff --git a/examples/CoolingBox/cooling_box_test.c b/examples/CoolingBox/cooling_box_test.c
deleted file mode 100644
index af01c1246854dc771088666d74823caad2e4f407..0000000000000000000000000000000000000000
--- a/examples/CoolingBox/cooling_box_test.c
+++ /dev/null
@@ -1,213 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (C) 2015 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/>.
- *
- ******************************************************************************/
-
-#include <unistd.h>
-#include "cooling.h"
-#include "cooling_struct.h"
-#include "hydro.h"
-#include "physical_constants.h"
-#include "swift.h"
-#include "units.h"
-
-/*
- * @brief Assign particle density and entropy corresponding to the
- * hydrogen number density and internal energy specified.
- *
- * @param p Particle data structure
- * @param cooling Cooling function data structure
- * @param cosmo Cosmology data structure
- * @param phys_const Physical constants data structure
- * @param nh Hydrogen number density (cgs units)
- * @param u Internal energy (cgs units)
- */
-void set_quantities(struct part *restrict p,
- struct xpart *restrict xp,
- const struct unit_system *restrict us,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
- const struct phys_const *restrict phys_const, float nh_cgs,
- double u) {
-
- float scale_factor = 1.0 / (1.0 + cosmo->z);
- double hydrogen_number_density =
- nh_cgs * units_cgs_conversion_factor(us, UNIT_CONV_LENGTH) *
- units_cgs_conversion_factor(us, UNIT_CONV_LENGTH) *
- units_cgs_conversion_factor(us, UNIT_CONV_LENGTH);
- p->rho = hydrogen_number_density * phys_const->const_proton_mass /
- p->chemistry_data.metal_mass_fraction[chemistry_element_H] * (cosmo->a * cosmo->a * cosmo->a);
-
- float pressure = (u * scale_factor * scale_factor) / cooling->internal_energy_scale *
- p->rho * (hydro_gamma_minus_one);
- p->entropy = pressure * (pow(p->rho, -hydro_gamma));
- xp->entropy_full = p->entropy;
-
- // Using hydro_set_init_internal_energy seems to work better for higher z for
- // setting the internal energy correctly However, with Gadget2 this just sets
- // the entropy to the internal energy, which needs to be converted somehow
- if (cosmo->z >= 1)
- hydro_set_init_internal_energy(
- p, (u * scale_factor * scale_factor) / cooling->internal_energy_scale);
-}
-
-/*
- * @brief Produces contributions to cooling rates for different
- * hydrogen number densities, from different metals,
- * tests 1d and 4d table interpolations produce
- * same results for cooling rate, dlambda/du and temperature.
- */
-int main(int argc, char **argv) {
- // Declare relevant structs
- struct swift_params *params = malloc(sizeof(struct swift_params));
- struct unit_system us;
- struct chemistry_global_data chem_data;
- struct part p;
- struct xpart xp;
- struct phys_const phys_const;
- struct cooling_function_data cooling;
- struct cosmology cosmo;
- char *parametersFileName = "./coolingBox.yml";
-
- float nh; // hydrogen number density
- double u; // internal energy
- const int npts = 250; // number of values for the internal energy at which
- // cooling rate is evaluated
-
- // Read options
- int param;
- float redshift = -1.0,
- log_10_nh =
- 100; // unreasonable values will be changed if not set in options
- while ((param = getopt(argc, argv, "z:d:m:t")) != -1) switch (param) {
- case 'z':
- // read redshift
- redshift = atof(optarg);
- break;
- case 'd':
- // read log10 of hydrogen number density
- log_10_nh = atof(optarg);
- break;
- case 'm':
- // read which yml file we need to use
- parametersFileName = optarg;
- break;
- case '?':
- if (optopt == 'z')
- printf("Option -%c requires an argument.\n", optopt);
- else
- printf("Unknown option character `\\x%x'.\n", optopt);
- error("invalid option(s) to testCooling");
- }
-
- // Read the parameter file
- if (params == NULL) error("Error allocating memory for the parameter file.");
- message("Reading runtime parameters from file '%s'", parametersFileName);
- parser_read_file(parametersFileName, params);
-
- // Init units
- units_init_from_params(&us, params, "InternalUnitSystem");
- phys_const_init(&us, params, &phys_const);
-
- // Init chemistry
- chemistry_init(params, &us, &phys_const, &chem_data);
- chemistry_first_init_part(&phys_const, &us, &cosmo, &chem_data, &p, &xp);
- chemistry_print(&chem_data);
-
- // Init cosmology
- cosmology_init(params, &us, &phys_const, &cosmo);
- cosmology_print(&cosmo);
- if (redshift == -1.0) {
- cosmo.z = 7.0;
- } else {
- cosmo.z = redshift;
- }
- message("redshift %.5e", cosmo.z);
-
- // Init cooling
- cooling_init(params, &us, &phys_const, &cooling);
- cooling_print(&cooling);
- cooling_update(&cosmo, &cooling, 0);
-
- // Calculate abundance ratios
- float *abundance_ratio;
- abundance_ratio = malloc((chemistry_element_count + 2) * sizeof(float));
- abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
-
- // extract mass fractions, calculate table indices and offsets
- float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
- float HeFrac =
- p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
- (XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
- int He_i, n_h_i;
- float d_He, d_n_h;
- get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
-
- // Calculate contributions from metals to cooling rate
- // open file
- const char *output_filename = "cooling_output.dat";
- FILE *output_file = fopen(output_filename, "w");
- if (output_file == NULL) {
- printf("Error opening file!\n");
- exit(1);
- }
-
- // set hydrogen number density
- if (log_10_nh == 100) {
- // hydrogen number density not specified in options
- nh = 1.0e-1;
- } else {
- nh = exp(M_LN10 * log_10_nh);
- }
-
- // set internal energy to dummy value, will get reset when looping over
- // internal energies
- u = 1.0e14;
- set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh, u);
- float inn_h = hydro_get_physical_density(&p, &cosmo) * XH / phys_const.const_proton_mass
- *cooling.number_density_scale;
- get_index_1d(cooling.nH, cooling.N_nH, log10(inn_h), &n_h_i, &d_n_h);
- message("inn_h %.5e nh %.5e", inn_h, nh);
-
- // Loop over internal energy
- float du;
- double temperature;
- for (int j = 0; j < npts; j++) {
- set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh,
- exp(M_LN10 * (10.0 + j * 8.0 / npts)));
- u = hydro_get_physical_internal_energy(&p, &xp, &cosmo) *
- cooling.internal_energy_scale;
- float cooling_du_dt;
-
- // calculate cooling rates
- double dLambdaNet_du;
- temperature = eagle_convert_u_to_temp(log10(u), &du, n_h_i, He_i, d_n_h, d_He, &cooling, &cosmo);
- //cooling_du_dt = eagle_print_metal_cooling_rate(
- // n_h_i, d_n_h, He_i, d_He, &p, &xp, &cooling, &cosmo, &phys_const,
- // abundance_ratio);
- cooling_du_dt = eagle_cooling_rate(
- log(u), &dLambdaNet_du, n_h_i, d_n_h, He_i, d_He, &p, &cooling, &cosmo, &phys_const,
- abundance_ratio);
- fprintf(output_file, "%.5e %.5e\n", u, cooling_du_dt);
- }
- fclose(output_file);
- message("done cooling rates test");
-
- free(params);
- return 0;
-}
-
diff --git a/examples/CoolingBox/cooling_rates_plot.py b/examples/CoolingBox/cooling_rates_plot.py
deleted file mode 100644
index 663cab522de9f5d11337347c360c947c232ab10b..0000000000000000000000000000000000000000
--- a/examples/CoolingBox/cooling_rates_plot.py
+++ /dev/null
@@ -1,129 +0,0 @@
-# Plots contribution to cooling rates from each of the different metals
-# based on cooling_output.dat and cooling_element_*.dat files produced
-# by testCooling.
-
-import matplotlib.pyplot as plt
-import numpy as np
-
-k_in_J_K = 1.38064852e-23
-mH_in_kg = 1.6737236e-27
-erg_to_J = 1.0e-7
-gas_gamma=5.0/3.0
-
-# Primoridal ean molecular weight as a function of temperature
-def mu(T):
- H_frac=0.752
- T_trans=10000.0
- 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=0.752
- T_trans=10000.0
- 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))
- if np.sum(mask_ionized) > 0:
- ret[mask_ionized] = T_over_mu[mask_ionized] * mu(T_trans*10)
-
- # Neutral gas?
- mask_neutral = (T_over_mu < (T_trans-1) / mu((T_trans-1)))
- if np.sum(mask_neutral) > 0:
- ret[mask_neutral] = T_over_mu[mask_neutral] * mu(0)
-
- return ret
-
-# Number of metals tracked by EAGLE cooling
-elements = 11
-
-# Declare arrays of internal energy and cooling rate
-u = [[] for i in range(elements+1)]
-cooling_rate = [[] for i in range(elements+1)]
-Temperature = [[] for i in range(elements+1)]
-
-# Read in total cooling rate
-#file_in = open('cooling_output.dat', 'r')
-#for line in file_in:
-# data = line.split()
-# u.append(float(data[0]))
-# cooling_rate[0].append(-float(data[1]))
-#
-#file_in.close()
-
-n_iz = 10
-for iz in range(n_iz):
- file_in = open('cooling_element_'+str(iz)+'.dat','r')
- z = float(file_in.readline())
- for line in file_in:
- data = line.split()
- u[iz+1].append(float(data[0]))
- cooling_rate[iz+1].append(-float(data[1]))
- file_in.close()
- a = 1.0/(1.0 + z)
- #u[iz+1] = np.asarray(u[iz+1]) / a**(3 * (gas_gamma - 1.))
- #Temperature[iz+1] = T(u[iz+1] * erg_to_J)
-
-# Plot
-ax = plt.subplot(111)
-#p0, = plt.loglog(u, cooling_rate[0], linewidth = 0.5, color = 'k', label = 'Total')
-p1, = plt.loglog(u[1], cooling_rate[1], linewidth = 0.5, color = 'k', label = 'z = 10')
-p2, = plt.loglog(u[2], cooling_rate[2], linewidth = 0.5, color = 'b', label = 'z = 9')
-p3, = plt.loglog(u[3], cooling_rate[3], linewidth = 0.5, color = 'g', label = 'z = 8')
-p4, = plt.loglog(u[4], cooling_rate[4], linewidth = 0.5, color = 'r', label = 'z = 7')
-p5, = plt.loglog(u[5], cooling_rate[5], linewidth = 0.5, color = 'c', label = 'z = 6')
-p6, = plt.loglog(u[6], cooling_rate[6], linewidth = 0.5, color = 'm', label = 'z = 5')
-p7, = plt.loglog(u[7], cooling_rate[7], linewidth = 0.5, color = 'y', label = 'z = 4')
-p8, = plt.loglog(u[8], cooling_rate[8], linewidth = 0.5, color = 'lightgray', label = 'z = 3')
-p9, = plt.loglog(u[9], cooling_rate[9], linewidth = 0.5, color = 'olive', label = 'z = 2')
-p10, = plt.loglog(u[10], cooling_rate[10], linewidth = 0.5, color = 'saddlebrown', label = 'z = 1')
-p11, = plt.loglog(u[1], -np.asarray(cooling_rate[1]), linewidth = 0.5, linestyle = '--', color = 'k')
-p12, = plt.loglog(u[2], -np.asarray(cooling_rate[2]), linewidth = 0.5, linestyle = '--', color = 'b')
-p13, = plt.loglog(u[3], -np.asarray(cooling_rate[3]), linewidth = 0.5, linestyle = '--', color = 'g')
-p14, = plt.loglog(u[4], -np.asarray(cooling_rate[4]), linewidth = 0.5, linestyle = '--', color = 'r')
-p15, = plt.loglog(u[5], -np.asarray(cooling_rate[5]), linewidth = 0.5, linestyle = '--', color = 'c')
-p16, = plt.loglog(u[6], -np.asarray(cooling_rate[6]), linewidth = 0.5, linestyle = '--', color = 'm')
-p17, = plt.loglog(u[7], -np.asarray(cooling_rate[7]), linewidth = 0.5, linestyle = '--', color = 'y')
-p18, = plt.loglog(u[8], -np.asarray(cooling_rate[8]), linewidth = 0.5, linestyle = '--', color = 'lightgray')
-p19, = plt.loglog(u[9], -np.asarray(cooling_rate[9]), linewidth = 0.5, linestyle = '--', color = 'olive')
-p20, = plt.loglog(u[10], -np.asarray(cooling_rate[10]), linewidth = 0.5, linestyle = '--', color = 'saddlebrown')
-ax.set_position([0.15,0.15,0.75,0.75])
-#plt.xlim([1e10,1e17])
-#plt.ylim([1e-24,1e-21])
-plt.xlabel("Internal energy ${\\rm{[erg \cdot g^{-1}]}}$", fontsize = 14)
-plt.ylabel("${\Lambda/n_h^2 }$ ${\\rm{[erg \cdot cm^3 \cdot s^{-1}]}}$", fontsize = 14)
-plt.legend(handles = [p1,p2,p3,p4,p5,p6,p7,p8,p9,p10])
-plt.show()
-
-# Read in contributions to cooling rates from each of the elements
-#for elem in range(elements):
-# file_in = open('cooling_element_'+str(elem)+'.dat','r')
-# for line in file_in:
-# data = line.split()
-# cooling_rate[elem+1].append(-float(data[0]))
-# file_in.close()
-
-## Plot
-#ax = plt.subplot(111)
-#p0, = plt.loglog(u, cooling_rate[0], linewidth = 0.5, color = 'k', label = 'Total')
-#p1, = plt.loglog(u, cooling_rate[1], linewidth = 0.5, color = 'k', linestyle = '--', label = 'H + He')
-#p2, = plt.loglog(u, cooling_rate[3], linewidth = 0.5, color = 'b', label = 'Carbon')
-#p3, = plt.loglog(u, cooling_rate[4], linewidth = 0.5, color = 'g', label = 'Nitrogen')
-#p4, = plt.loglog(u, cooling_rate[5], linewidth = 0.5, color = 'r', label = 'Oxygen')
-#p5, = plt.loglog(u, cooling_rate[6], linewidth = 0.5, color = 'c', label = 'Neon')
-#p6, = plt.loglog(u, cooling_rate[7], linewidth = 0.5, color = 'm', label = 'Magnesium')
-#p7, = plt.loglog(u, cooling_rate[8], linewidth = 0.5, color = 'y', label = 'Silicon')
-#p8, = plt.loglog(u, cooling_rate[9], linewidth = 0.5, color = 'lightgray', label = 'Sulphur')
-#p9, = plt.loglog(u, cooling_rate[10], linewidth = 0.5, color = 'olive', label = 'Calcium')
-#p10, = plt.loglog(u, cooling_rate[11], linewidth = 0.5, color = 'saddlebrown', label = 'Iron')
-#ax.set_position([0.15,0.15,0.75,0.75])
-#plt.xlim([1e12,1e17])
-#plt.ylim([1e-24,1e-21])
-#plt.xlabel("Internal energy ${\\rm{[erg \cdot g^{-1}]}}$", fontsize = 14)
-#plt.ylabel("${\Lambda/n_h^2 }$ ${\\rm{[erg \cdot cm^3 \cdot s^{-1}]}}$", fontsize = 14)
-##plt.legend(handles = [p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10])
-#plt.show()
diff --git a/examples/CoolingBox/getGlass.sh b/examples/CoolingBox/getGlass.sh
index 81c2cc6ff2b3aa07cac3a6ee7a225f9e9a6b740a..ffd92e88deae6e91237059adac2a6c2067caee46 100755
--- a/examples/CoolingBox/getGlass.sh
+++ b/examples/CoolingBox/getGlass.sh
@@ -1,2 +1,2 @@
#!/bin/bash
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_16.hdf5
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5
diff --git a/examples/CoolingBox/makeIC.py b/examples/CoolingBox/makeIC.py
index a8e888526fb1f86dd1341af377c91a480f124ca5..f863e174b1fcd404ae178fe324c7a165598b4af0 100644
--- a/examples/CoolingBox/makeIC.py
+++ b/examples/CoolingBox/makeIC.py
@@ -26,21 +26,17 @@ from numpy import *
# Parameters
periodic= 1 # 1 For periodic box
-boxSize = 3.0857e21 # 1 kiloparsec
-rho = 2.36748e-25 # Density in cgs
-P = 6.68e-12 # Pressure in cgs
+boxSize = 1 # 1 kiloparsec
+rho = 3.2e3 # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
+P = 4.5e6 # Pressure in code units (at 10^5K)
gamma = 5./3. # Gas adiabatic index
eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "coolingBox.hdf5"
-if len(sys.argv) == 2:
- rho = float(sys.argv[1])
- P = float(sys.argv[2])
-
#---------------------------------------------------
# Read id, position and h from glass
-glass = h5py.File("glassCube_16.hdf5", "r")
+glass = h5py.File("glassCube_32.hdf5", "r")
ids = glass["/PartType0/ParticleIDs"][:]
pos = glass["/PartType0/Coordinates"][:,:] * boxSize
h = glass["/PartType0/SmoothingLength"][:] * boxSize
@@ -70,9 +66,9 @@ grp.attrs["PeriodicBoundariesOn"] = periodic
# Units
grp = file.create_group("/Units")
-grp.attrs["Unit length in cgs (U_L)"] = 1.
-grp.attrs["Unit mass in cgs (U_M)"] = 1.
-grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit length in cgs (U_L)"] = 3.0857e21
+grp.attrs["Unit mass in cgs (U_M)"] = 2.0e33
+grp.attrs["Unit time in cgs (U_t)"] = 3.0857e16
grp.attrs["Unit current in cgs (U_I)"] = 1.
grp.attrs["Unit temperature in cgs (U_T)"] = 1.
diff --git a/examples/CoolingBox/test.sh b/examples/CoolingBox/test.sh
deleted file mode 100755
index 8f746d22f194fc40502c24d693e50a75e3dd1ee1..0000000000000000000000000000000000000000
--- a/examples/CoolingBox/test.sh
+++ /dev/null
@@ -1,57 +0,0 @@
-#!/bin/bash
-
-# loop over redshift
-for z in 0.1; do
- # loop over solar and zero metal abundances
- for solar in 1; do
- # loop over log_10 hydrogen number density
- for nh_exp in -1; do
- # Choose which yml file to use for different metallicities
- if [ $solar == 1 ]
- then
- yml_file="coolingBox_solar.yml"
- elif [ $solar == 0 ]
- then
- yml_file="coolingBox_primordial.yml"
- fi
-
- # Calculate cooling rates for given redshift, density and metallicity
- rm cooling_*.dat
- ./testCooling -z $z -d $nh_exp -m $yml_file
-
-
- # set starting, ending redshift, how often to write to file
- a_begin=$(python -c "print 1.0/(1.0+$z)")
- a_end=$(python -c "print min(1.0, $a_begin*1.0001)")
- delta_a=$(python -c "print 1.0 + ($a_end/$a_begin - 1.0)/50.")
-
- # change hydrogen number density
- if [ $solar == 0 ];
- then
- rho=$(python -c "print 10.0**$nh_exp/0.7*1.6726e-24")
- else
- rho=$(python -c "print 10.0**$nh_exp/0.752*1.6726e-24")
- fi
- for pressure_index in 8; do
- pressure=$(python -c "print 6.68e-14")
-
- python makeIC.py $rho $pressure
- rm coolingBox_*hdf5
-
- # run cooling box
- ../swift -s -c -C -t 4 -P Cosmology:a_begin:$a_begin -P Cosmology:a_end:$a_end -P Snapshots:scale_factor_first:$a_begin -P Snapshots:delta_time:$delta_a -P Statistics:scale_factor_first:$a_begin -P Statistics:delta_time:$delta_a $yml_file 2>&1 | tee output.log
-
- max=0
- for file in $( ls -lth coolingBox_*.hdf5 | head -n 1 ); do
- f=${file//hdf5/}
- f=${f//[!0-9]/}
- done
- frames=$((10#$f))
- echo "number of frames $frames"
-
- # check if everything worked and create plots
- python analytical_test.py $z $nh_exp $pressure_index $solar $frames
- done
- done
- done
-done
diff --git a/examples/CoolingBox/testCooling.c b/examples/CoolingBox/testCooling.c
deleted file mode 100644
index 37db87f2bcbbb493128c3b7c7923c70502c4ed63..0000000000000000000000000000000000000000
--- a/examples/CoolingBox/testCooling.c
+++ /dev/null
@@ -1,187 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (C) 2015 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/>.
- *
- ******************************************************************************/
-
-#include <unistd.h>
-#include "cooling.h"
-#include "cooling_struct.h"
-#include "hydro.h"
-#include "physical_constants.h"
-#include "swift.h"
-#include "units.h"
-
-/*
- * @brief Assign particle density and entropy corresponding to the
- * hydrogen number density and internal energy specified.
- *
- * @param p Particle data structure
- * @param xp extra particle structure
- * @param us unit system struct
- * @param cooling Cooling function data structure
- * @param cosmo Cosmology data structure
- * @param phys_const Physical constants data structure
- * @param nh_cgs Hydrogen number density (cgs units)
- * @param u Internal energy (cgs units)
- * @param ti_current integertime to set cosmo quantities
- */
-void set_quantities(struct part *restrict p,
- struct xpart *restrict xp,
- const struct unit_system *restrict us,
- const struct cooling_function_data *restrict cooling,
- struct cosmology *restrict cosmo,
- const struct phys_const *restrict phys_const, float nh_cgs,
- double u,
- integertime_t ti_current) {
-
- /* Update cosmology quantities */
- cosmology_update(cosmo,phys_const,ti_current);
-
- /* calculate density */
- double hydrogen_number_density = nh_cgs / cooling->number_density_scale;
- p->rho = hydrogen_number_density * phys_const->const_proton_mass /
- p->chemistry_data.metal_mass_fraction[chemistry_element_H] * (cosmo->a * cosmo->a * cosmo->a);
-
- /* update entropy based on internal energy */
- float pressure = (u * cosmo->a * cosmo->a) / cooling->internal_energy_scale *
- p->rho * (hydro_gamma_minus_one);
- p->entropy = pressure * (pow(p->rho, -hydro_gamma));
- xp->entropy_full = p->entropy;
-
-}
-
-/*
- * @brief Produces contributions to cooling rates for different
- * hydrogen number densities, from different metals,
- * tests 1d and 4d table interpolations produce
- * same results for cooling rate, dlambda/du and temperature.
- */
-int main(int argc, char **argv) {
- // Declare relevant structs
- struct swift_params *params = malloc(sizeof(struct swift_params));
- struct unit_system us;
- struct chemistry_global_data chem_data;
- struct part p;
- struct xpart xp;
- struct phys_const phys_const;
- struct cooling_function_data cooling;
- struct cosmology cosmo;
- char *parametersFileName = "./coolingBox_solar.yml";
-
- float nh; // hydrogen number density
- double u; // internal energy
-
- /* Number of values to test for in redshift,
- * hydrogen number density and internal energy */
- const int n_z = 50;
- const int n_nh = 50;
- const int n_u = 50;
-
- /* Number of subcycles and tolerance used to compare
- * subcycled and implicit solution. Note, high value
- * of tolerance due to mismatch between explicit and
- * implicit solution for large timesteps */
- const int n_subcycle = 1000;
- const float integration_tolerance = 0.2;
-
- /* Set dt */
- const float dt_cool = 1.0e-5;
- const float dt_therm = 1.0e-5;
-
- /* Read the parameter file */
- if (params == NULL) error("Error allocating memory for the parameter file.");
- message("Reading runtime parameters from file '%s'", parametersFileName);
- parser_read_file(parametersFileName, params);
-
- /* Init units */
- units_init_from_params(&us, params, "InternalUnitSystem");
- phys_const_init(&us, params, &phys_const);
-
- /* Init chemistry */
- chemistry_init(params, &us, &phys_const, &chem_data);
- chemistry_first_init_part(&phys_const, &us, &cosmo, &chem_data, &p, &xp);
- chemistry_print(&chem_data);
-
- /* Init cosmology */
- cosmology_init(params, &us, &phys_const, &cosmo);
- cosmology_print(&cosmo);
-
- /* Init cooling */
- cooling_init(params, &us, &phys_const, &cooling);
- cooling_print(&cooling);
- cooling_update(&cosmo, &cooling, 0);
-
- /* Calculate abundance ratios */
- float *abundance_ratio;
- abundance_ratio = malloc((chemistry_element_count + 2) * sizeof(float));
- abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
-
- /* extract mass fractions, calculate table indices and offsets */
- float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
- float HeFrac =
- p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
- (XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
- int He_i;
- float d_He;
- get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
-
- /* Cooling function needs to know the minimal energy. Set it to the lowest
- * internal energy in the cooling table. */
- struct hydro_props hydro_properties;
- hydro_properties.minimal_internal_energy = exp(M_LN10 * cooling.Therm[0]) / cooling.internal_energy_scale;
-
- /* calculate spacing in nh and u */
- const float delta_nh = (cooling.nH[cooling.N_nH - 1] - cooling.nH[0])/n_nh;
- const float delta_u = (cooling.Therm[cooling.N_Temp - 1] - cooling.Therm[0])/n_u;
-
- for(int z_i = 0; z_i < n_z; z_i++) {
- integertime_t ti_current = max_nr_timesteps/n_z*z_i;
- for(int nh_i = 0; nh_i < n_nh; nh_i++) {
- nh = exp(M_LN10 * cooling.nH[0] + delta_nh*nh_i);
- for (int u_i = 0; u_i < n_u; u_i++) {
- message("z_i %d n_h_i %d u_i %d", z_i, nh_i, u_i);
- u = exp(M_LN10 * cooling.Therm[0] + delta_u*u_i);
-
- /* update nh, u, z */
- set_quantities(&p,&xp,&us,&cooling,&cosmo,&phys_const, nh, u, ti_current);
-
- /* calculate subcycled solution */
- for (int t_subcycle = 0; t_subcycle < n_subcycle; t_subcycle++) {
- cooling_cool_part(&phys_const,&us,&cosmo,&hydro_properties,&cooling,&p,&xp,dt_cool/n_subcycle,dt_therm/n_subcycle);
- xp.entropy_full += p.entropy_dt*dt_therm/n_subcycle;
- }
- double u_subcycled = hydro_get_physical_internal_energy(&p,&xp,&cosmo)*cooling.internal_energy_scale;
-
- /* reset quantities to nh, u, and z that we want to test */
- set_quantities(&p,&xp,&us,&cooling,&cosmo,&phys_const, nh, u, ti_current);
-
- /* compute implicit solution */
- cooling_cool_part(&phys_const,&us,&cosmo,&hydro_properties,&cooling,&p,&xp,dt_cool,dt_therm);
- double u_implicit = hydro_get_physical_internal_energy(&p,&xp,&cosmo)*cooling.internal_energy_scale;
-
- /* check if the two solutions are consistent */
- if (fabs((u_implicit-u_subcycled)/u_subcycled) > integration_tolerance) message("implicit and subcycled solutions do not match. z_i %d nh_i %d u_i %d implicit %.5e subcycled %.5e error %.5e", z_i, nh_i, u_i, u_implicit, u_subcycled, fabs((u_implicit-u_subcycled)/u_subcycled));
- }
- }
- }
-
- message("done test");
-
- free(params);
- return 0;
-}
-
diff --git a/examples/CoolingBox/testCooling_cool_part.c b/examples/CoolingBox/testCooling_cool_part.c
deleted file mode 100644
index c2529e1f196edc6803f52890a8df0886c301bf08..0000000000000000000000000000000000000000
--- a/examples/CoolingBox/testCooling_cool_part.c
+++ /dev/null
@@ -1,190 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (C) 2015 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/>.
- *
- ******************************************************************************/
-
-#include <unistd.h>
-#include "cooling.h"
-#include "cooling_struct.h"
-#include "hydro.h"
-#include "physical_constants.h"
-#include "swift.h"
-#include "units.h"
-
-/*
- * @brief Assign particle density and entropy corresponding to the
- * hydrogen number density and internal energy specified.
- *
- * @param p Particle data structure
- * @param cooling Cooling function data structure
- * @param cosmo Cosmology data structure
- * @param phys_const Physical constants data structure
- * @param nh Hydrogen number density (cgs units)
- * @param u Internal energy (cgs units)
- */
-void set_quantities(struct part *restrict p,
- struct xpart *restrict xp,
- const struct unit_system *restrict us,
- const struct cooling_function_data *restrict cooling,
- struct cosmology *restrict cosmo,
- const struct phys_const *restrict phys_const, float nh_cgs,
- double u,
- float z) {
-
- cosmo->z = z;
- float scale_factor = 1.0 / (1.0 + cosmo->z);
- double hydrogen_number_density =
- nh_cgs * units_cgs_conversion_factor(us, UNIT_CONV_LENGTH) *
- units_cgs_conversion_factor(us, UNIT_CONV_LENGTH) *
- units_cgs_conversion_factor(us, UNIT_CONV_LENGTH);
- p->rho = hydrogen_number_density * phys_const->const_proton_mass /
- p->chemistry_data.metal_mass_fraction[chemistry_element_H] * (scale_factor * scale_factor * scale_factor);
-
- float pressure = (u * scale_factor * scale_factor) / cooling->internal_energy_scale *
- p->rho * (hydro_gamma_minus_one);
- p->entropy = pressure * (pow(p->rho, -hydro_gamma));
- xp->entropy_full = p->entropy;
-
- // Using hydro_set_init_internal_energy seems to work better for higher z for
- // setting the internal energy correctly However, with Gadget2 this just sets
- // the entropy to the internal energy, which needs to be converted somehow
- if (cosmo->z >= 1)
- hydro_set_init_internal_energy(
- p, (u * scale_factor * scale_factor) / cooling->internal_energy_scale);
-}
-
-/*
- * @brief Produces contributions to cooling rates for different
- * hydrogen number densities, from different metals,
- * tests 1d and 4d table interpolations produce
- * same results for cooling rate, dlambda/du and temperature.
- */
-int main(int argc, char **argv) {
- // Declare relevant structs
- struct swift_params *params = malloc(sizeof(struct swift_params));
- struct unit_system us;
- struct chemistry_global_data chem_data;
- struct part p;
- struct xpart xp;
- struct phys_const phys_const;
- struct cooling_function_data cooling;
- struct cosmology cosmo;
- char *parametersFileName = "./coolingBox.yml";
-
- float nh; // hydrogen number density
- double u; // internal energy
-
- // Read options
- int param;
- float redshift = -1.0,
- log_10_nh =
- 100; // unreasonable values will be changed if not set in options
- while ((param = getopt(argc, argv, "z:d:m:t")) != -1) switch (param) {
- case 'z':
- // read redshift
- redshift = atof(optarg);
- break;
- case 'd':
- // read log10 of hydrogen number density
- log_10_nh = atof(optarg);
- break;
- case 'm':
- // read which yml file we need to use
- parametersFileName = optarg;
- break;
- case '?':
- if (optopt == 'z')
- printf("Option -%c requires an argument.\n", optopt);
- else
- printf("Unknown option character `\\x%x'.\n", optopt);
- error("invalid option(s) to testCooling");
- }
-
- // Read the parameter file
- if (params == NULL) error("Error allocating memory for the parameter file.");
- message("Reading runtime parameters from file '%s'", parametersFileName);
- parser_read_file(parametersFileName, params);
-
- // Init units
- units_init_from_params(&us, params, "InternalUnitSystem");
- phys_const_init(&us, params, &phys_const);
-
- // Init chemistry
- chemistry_init(params, &us, &phys_const, &chem_data);
- chemistry_first_init_part(&phys_const, &us, &cosmo, &chem_data, &p, &xp);
- chemistry_print(&chem_data);
-
- // Init cosmology
- cosmology_init(params, &us, &phys_const, &cosmo);
- cosmology_print(&cosmo);
- if (redshift == -1.0) {
- cosmo.z = 7.0;
- } else {
- cosmo.z = redshift;
- }
- message("redshift %.5e", cosmo.z);
-
- // Init cooling
- cooling_init(params, &us, &phys_const, &cooling);
- cooling_print(&cooling);
- cooling_update(&cosmo, &cooling, 0);
-
- // Calculate abundance ratios
- float *abundance_ratio;
- abundance_ratio = malloc((chemistry_element_count + 2) * sizeof(float));
- abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
-
- // extract mass fractions, calculate table indices and offsets
- float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
- float HeFrac =
- p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
- (XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
- int He_i;//, n_h_i;
- float d_He;//, d_n_h;
- get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
-
- int n_z = 100;
- int n_nh = 100;
- int n_u = 100;
-
- struct hydro_props hydro_properties;
- hydro_properties.minimal_internal_energy = exp(M_LN10 * cooling.Therm[0]) / cooling.internal_energy_scale;
- float dt_cool = 0.1;
- float dt_therm = 0.1;
-
- float delta_z = (cooling.Redshifts[cooling.N_Redshifts - 1] - cooling.Redshifts[0])/n_z;
- float delta_nh = (cooling.nH[cooling.N_nH - 1] - cooling.nH[0])/n_nh;
- float delta_u = (cooling.Therm[cooling.N_Temp - 1] - cooling.Therm[0])/n_u;
- for(int i = 0; i < n_z; i++) {
- for(int j = 0; j < n_nh; j++) {
- for (int k = 0; k < n_u; k++) {
- float z = cooling.Redshifts[0] + delta_z*i;
- nh = exp(M_LN10 * cooling.nH[0] + delta_nh*j);
- u = exp(M_LN10 * cooling.Therm[0] + delta_u*k);
- set_quantities(&p,&xp,&us,&cooling,&cosmo,&phys_const, nh, u, z);
-
- cooling_cool_part(&phys_const,&us,&cosmo,&hydro_properties,&cooling,&p,&xp,dt_cool,dt_therm);
- }
- }
- }
-
- message("done test");
-
- free(params);
- return 0;
-}
-
diff --git a/examples/CoolingBox/test_script.sh b/examples/CoolingBox/test_script.sh
deleted file mode 100644
index 119af147ece06c365095290f97647734f9dbd13b..0000000000000000000000000000000000000000
--- a/examples/CoolingBox/test_script.sh
+++ /dev/null
@@ -1,116 +0,0 @@
-#!/bin/bash
-
-max_number() {
- printf "%s\n" "$@" | sort -g | tail -n1
-}
-
-# loop over redshift
-for z in 6.0; do
- # loop over solar and zero metal abundances
- for solar in 0; do
- # loop over log_10 hydrogen number density
- for nh_exp in $(seq -4 1 -1); do
- #change parameters in yml file for calculating explicit ode solution of cooling
- cd ../CoolingRates
- if [ $solar == 1 ]
- then
- sed -i "/InitMetallicity: / s/: \+[[:alnum:],\.,-]\+/: 0.014/g" testCooling.yml
- sed -i "/InitAbundance_Hydrogen: / s/: \+[[:alnum:],\.,-]\+/: 0.70649785/g" testCooling.yml
- sed -i "/InitAbundance_Helium: / s/: \+[[:alnum:],\.,-]\+/: 0.28055534/g" testCooling.yml
- sed -i "/InitAbundance_Carbon: / s/: \+[[:alnum:],\.,-]\+/: 2.0665436e-3/g" testCooling.yml
- sed -i "/InitAbundance_Nitrogen: / s/: \+[[:alnum:],\.,-]\+/: 8.3562563e-4/g" testCooling.yml
- sed -i "/InitAbundance_Oxygen: / s/: \+[[:alnum:],\.,-]\+/: 5.4926244e-3/g" testCooling.yml
- sed -i "/InitAbundance_Neon: / s/: \+[[:alnum:],\.,-]\+/: 1.4144605e-3/g" testCooling.yml
- sed -i "/InitAbundance_Magnesium: / s/: \+[[:alnum:],\.,-]\+/: 5.907064e-4/g" testCooling.yml
- sed -i "/InitAbundance_Silicon: / s/: \+[[:alnum:],\.,-]\+/: 6.825874e-4/g" testCooling.yml
- sed -i "/InitAbundance_Iron: / s/: \+[[:alnum:],\.,-]\+/: 1.1032152e-3/g" testCooling.yml
- sed -i "/CalciumOverSilicon: / s/: \+[[:alnum:],\.,-]\+/: 0.0941736/g" testCooling.yml
- sed -i "/SulphurOverSilicon: / s/: \+[[:alnum:],\.,-]\+/: 0.6054160/g" testCooling.yml
- elif [ $solar == 0 ]
- then
- sed -i "/InitMetallicity: / s/: \+[[:alnum:],\.,-]\+/: 0.0/g" testCooling.yml
- sed -i "/InitAbundance_Hydrogen: / s/: \+[[:alnum:],\.,-]\+/: 0.752/g" testCooling.yml
- sed -i "/InitAbundance_Helium: / s/: \+[[:alnum:],\.,-]\+/: 0.248/g" testCooling.yml
- sed -i "/InitAbundance_Carbon: / s/: \+[[:alnum:],\.,-]\+/: 0.000/g" testCooling.yml
- sed -i "/InitAbundance_Nitrogen: / s/: \+[[:alnum:],\.,-]\+/: 0.000/g" testCooling.yml
- sed -i "/InitAbundance_Oxygen: / s/: \+[[:alnum:],\.,-]\+/: 0.000/g" testCooling.yml
- sed -i "/InitAbundance_Neon: / s/: \+[[:alnum:],\.,-]\+/: 0.000/g" testCooling.yml
- sed -i "/InitAbundance_Magnesium: / s/: \+[[:alnum:],\.,-]\+/: 0.000/g" testCooling.yml
- sed -i "/InitAbundance_Silicon: / s/: \+[[:alnum:],\.,-]\+/: 0.000/g" testCooling.yml
- sed -i "/InitAbundance_Iron: / s/: \+[[:alnum:],\.,-]\+/: 0.000/g" testCooling.yml
- sed -i "/CalciumOverSilicon: / s/: \+[[:alnum:],\.,-]\+/: 0.0941736/g" testCooling.yml
- sed -i "/SulphurOverSilicon: / s/: \+[[:alnum:],\.,-]\+/: 0.6054160/g" testCooling.yml
- fi
- rm cooling_*.dat
- ./testCooling -z $z -d $nh_exp
- cd ../CoolingBox
- cp ../CoolingRates/cooling_output.dat ./
-
- #change parameters in coolingBox.yml
- if [ $solar == 1 ]
- then
- sed -i "/InitMetallicity: / s/: \+[[:alnum:],\.,-]\+/: 0.014/g" coolingBox.yml
- sed -i "/InitAbundance_Hydrogen: / s/: \+[[:alnum:],\.,-]\+/: 0.70649785/g" coolingBox.yml
- sed -i "/InitAbundance_Helium: / s/: \+[[:alnum:],\.,-]\+/: 0.28055534/g" coolingBox.yml
- sed -i "/InitAbundance_Carbon: / s/: \+[[:alnum:],\.,-]\+/: 2.0665436e-3/g" coolingBox.yml
- sed -i "/InitAbundance_Nitrogen: / s/: \+[[:alnum:],\.,-]\+/: 8.3562563e-4/g" coolingBox.yml
- sed -i "/InitAbundance_Oxygen: / s/: \+[[:alnum:],\.,-]\+/: 5.4926244e-3/g" coolingBox.yml
- sed -i "/InitAbundance_Neon: / s/: \+[[:alnum:],\.,-]\+/: 1.4144605e-3/g" coolingBox.yml
- sed -i "/InitAbundance_Magnesium: / s/: \+[[:alnum:],\.,-]\+/: 5.907064e-4/g" coolingBox.yml
- sed -i "/InitAbundance_Silicon: / s/: \+[[:alnum:],\.,-]\+/: 6.825874e-4/g" coolingBox.yml
- sed -i "/InitAbundance_Iron: / s/: \+[[:alnum:],\.,-]\+/: 1.1032152e-3/g" coolingBox.yml
- sed -i "/CalciumOverSilicon: / s/: \+[[:alnum:],\.,-]\+/: 0.0941736/g" coolingBox.yml
- sed -i "/SulphurOverSilicon: / s/: \+[[:alnum:],\.,-]\+/: 0.6054160/g" coolingBox.yml
- elif [ $solar == 0 ]
- then
- sed -i "/InitMetallicity: / s/: \+[[:alnum:],\.,-]\+/: 0.0/g" coolingBox.yml
- sed -i "/InitAbundance_Hydrogen: / s/: \+[[:alnum:],\.,-]\+/: 0.752/g" coolingBox.yml
- sed -i "/InitAbundance_Helium: / s/: \+[[:alnum:],\.,-]\+/: 0.248/g" coolingBox.yml
- sed -i "/InitAbundance_Carbon: / s/: \+[[:alnum:],\.,-]\+/: 0.000/g" coolingBox.yml
- sed -i "/InitAbundance_Nitrogen: / s/: \+[[:alnum:],\.,-]\+/: 0.000/g" coolingBox.yml
- sed -i "/InitAbundance_Oxygen: / s/: \+[[:alnum:],\.,-]\+/: 0.000/g" coolingBox.yml
- sed -i "/InitAbundance_Neon: / s/: \+[[:alnum:],\.,-]\+/: 0.000/g" coolingBox.yml
- sed -i "/InitAbundance_Magnesium: / s/: \+[[:alnum:],\.,-]\+/: 0.000/g" coolingBox.yml
- sed -i "/InitAbundance_Silicon: / s/: \+[[:alnum:],\.,-]\+/: 0.000/g" coolingBox.yml
- sed -i "/InitAbundance_Iron: / s/: \+[[:alnum:],\.,-]\+/: 0.000/g" coolingBox.yml
- sed -i "/CalciumOverSilicon: / s/: \+[[:alnum:],\.,-]\+/: 0.0941736/g" coolingBox.yml
- sed -i "/SulphurOverSilicon: / s/: \+[[:alnum:],\.,-]\+/: 0.6054160/g" coolingBox.yml
- fi
-
- # set starting, ending redshift, how often to write to file
- a_begin=$(python -c "print 1.0/(1.0+$z)")
- a_end=$(python -c "print min(1.0, $a_begin*1.01)")
- first_ouput_a=$a_begin
- delta_a=$(python -c "print 1.0 + ($a_end/$a_begin - 1.0)/500.")
- sed -i "/a_begin: / s/: \+[[:alnum:],\.,-]\+/: $a_begin/g" coolingBox.yml
- sed -i "/a_end: / s/: \+[[:alnum:],\.,-]\+/: $a_end/g" coolingBox.yml
- sed -i "/scale_factor_first: / s/: \+[[:alnum:],\.,-]\+/: $first_ouput_a/g" coolingBox.yml
- sed -i "/delta_time: / s/: \+[[:alnum:],\.,-]\+/: $delta_a/g" coolingBox.yml
-
- # change hydrogen number density
- nh=$(python -c "print 3.555*10.0**($nh_exp+7)/(1.0+$z)**3")
- sed -i "/^rho =/ s/= \S*/= $nh/g" makeIC.py
- for pressure_index in 6; do
- pressure=$(python -c "print 4.5*10.0**($pressure_index + $nh_exp + 3)/(1.0+$z)")
- # change pressure (and hence energy)
- sed -i "/^P =/ s/= \S*/= $pressure/g" makeIC.py
- python makeIC.py
- rm coolingBox_*hdf5
-
- # run cooling box
- ../swift -s -c -C -t 16 coolingBox.yml
-
- max=0
- for file in $( ls -lth coolingBox_*.hdf5 | head -n 1 ); do
- f=${file//hdf5/}
- f=${f//[!0-9]/}
- done
- frames=$((10#$f))
- echo "number of frames $frames"
-
- # check if everything worked and create plots
- python analytical_test.py $z $nh_exp $pressure_index $solar $frames
- done
- done
- done
-done
diff --git a/examples/SmoothedMetallicity/run.sh b/examples/SmoothedMetallicity/run.sh
index de8c55d678bcb611934af450940d8ed8e6c15d6b..590970d2d9c207df09c8b7890a6364996a321c1e 100755
--- a/examples/SmoothedMetallicity/run.sh
+++ b/examples/SmoothedMetallicity/run.sh
@@ -16,4 +16,4 @@ fi
../swift -n 1 -s -t 4 smoothed_metallicity.yml 2>&1 | tee output.log
# Plot the solution
-python plotSolution.py 1
+python plotSolution_eagle.py 1
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 9698e66d1a7a80566a9347bb629e9670eafbb98c..7a0242437b8496d8c8756b1bccd2abb4d991262f 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -41,7 +41,7 @@ check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC \
testGravityDerivatives testPotentialSelf testPotentialPair testEOS testUtilities \
testSelectOutput testCbrt testCosmology testOutputList test27cellsStars \
- test27cellsStars_subset
+ test27cellsStars_subset testCooling
# Rebuild tests when SWIFT is updated.
$(check_PROGRAMS): ../src/.libs/libswiftsim.a
@@ -128,6 +128,8 @@ testEOS_SOURCES = testEOS.c
testUtilities_SOURCES = testUtilities.c
+testCooling_SOURCES = testCooling.c
+
# Files necessary for distribution
EXTRA_DIST = testReading.sh makeInput.py testActivePair.sh \
test27cells.sh test27cellsPerturbed.sh testParser.sh testPeriodicBC.sh \