diff --git a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/advect_ions.yml b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/advect_ions.yml
index 4558c2ecd474c657cc069d8b74c5a4ac1ae48e96..89271db00a87eddb8f76f496b7e0e2486c9acd28 100644
--- a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/advect_ions.yml
+++ b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/advect_ions.yml
@@ -13,7 +13,7 @@ InternalUnitSystem:
TimeIntegration:
max_nr_rt_subcycles: 1
time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 0.01 # end time: radiation reaches edge of box
+ time_end: 1 # end time: radiation reaches edge of box
dt_min: 1.e-12 # The minimal time-step size of the simulation (in internal units).
dt_max: 1.e-03 # The maximal time-step size of the simulation (in internal units).
@@ -22,7 +22,7 @@ TimeIntegration:
Snapshots:
basename: output # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
- delta_time: 0.001 # Time between snapshots
+ delta_time: 0.1 # Time between snapshots
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/advect_ions_cosmo.yml b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/advect_ions_cosmo.yml
new file mode 100644
index 0000000000000000000000000000000000000000..e4debac1573ff0678b604b2670b5e007b45fc562
--- /dev/null
+++ b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/advect_ions_cosmo.yml
@@ -0,0 +1,63 @@
+MetaData:
+ run_name: advect_ions
+
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.98841586e+33 # 1 M_Sol
+ UnitLength_in_cgs: 3.08567758e21 # kpc in cm
+ UnitVelocity_in_cgs: 1.e5 # km/s
+ UnitCurrent_in_cgs: 1.
+ UnitTemp_in_cgs: 1. # K
+
+# Parameters governing the time integration
+TimeIntegration:
+ max_nr_rt_subcycles: 1
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 0.01 # end time: radiation reaches edge of box
+ dt_min: 1.e-12 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1.e-03 # The maximal time-step size of the simulation (in internal units).
+
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: cosmo_output # Common part of the name of output files
+ output_list_on: 1
+ output_list: snapshot_times_high_redshift
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ time_first: 0.00990099
+ delta_time: 1.06 # 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.9 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 10. # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./advect_ions.hdf5 # The file to read
+ periodic: 1 # periodic ICs. Keep them periodic so we don't loose photon energy.
+
+GEARRT:
+ f_reduce_c: 1.e-5 # reduce the speed of light for the RT solver by multiplying c with this factor
+ CFL_condition: 0.9
+ f_limit_cooling_time: 0.0 # (Optional) multiply the cooling time by this factor when estimating maximal next time step. Set to 0.0 to turn computation of cooling time off.
+ photon_groups_Hz: [3.288e15, 5.945e15, 13.157e15] # Lower photon frequency group bin edges in Hz. Needs to have exactly N elements, where N is the configured number of bins --with-RT=GEAR_N
+ stellar_luminosity_model: const # Which model to use to determine the stellar luminosities.
+ const_stellar_luminosities_LSol: [0., 0., 0.] # (Conditional) constant star luminosities for each photon frequency group to use if stellar_luminosity_model:const is set, in units of Solar Luminosity.
+ hydrogen_mass_fraction: 0.76 # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
+ set_equilibrium_initial_ionization_mass_fractions: 0 # (Optional) set the initial ionization fractions depending on gas temperature assuming ionization equilibrium.
+ set_initial_ionization_mass_fractions: 0 # (Optional) manually overwrite initial mass fraction of each species (using the values you set below)
+ stellar_spectrum_type: 1 # Which radiation spectrum to use. 0: constant over all frequencies. 1: blackbody spectrum.
+ stellar_spectrum_blackbody_temperature_K: 1.e5 # (Conditional) if stellar_spectrum_type=1, use this temperature (in K) for the blackbody spectrum.
+ skip_thermochemistry: 1
+
+Cosmology: # Planck13 (EAGLE flavour)
+ a_begin: 0.00990099 # z~100
+ a_end: 0.01408451 # z~70
+ h: 0.6777
+ Omega_cdm: 0.2587481
+ Omega_lambda: 0.693
+ Omega_b: 0.0482519
diff --git a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/makeIC.py b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/makeIC.py
index 0f5621e826d5aa99ada36f38046e7fd70137d2ca..8c59ad795d0b360648fb9b35d47e285aa92f5457 100755
--- a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/makeIC.py
+++ b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/makeIC.py
@@ -43,7 +43,7 @@ if __name__ == "__main__":
# Set up metadata
unitL = unyt.Mpc
- edgelen = 2 * 15 * 1e-3 * unitL # 30 kpc
+ edgelen = 1 * unyt.Mpc
edgelen = edgelen.to(unitL)
boxsize = np.array([1.0, 1.0, 0.0]) * edgelen
diff --git a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/plotIonization.py b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/plotIonization.py
index 57f79fc99b2e333184ba342c531444fcf690574b..e76b9d6cbf1ebb85543ac09e7584a87380e4dec6 100755
--- a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/plotIonization.py
+++ b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/plotIonization.py
@@ -31,13 +31,12 @@ import swiftsimio
from matplotlib import pyplot as plt
from matplotlib.colors import SymLogNorm
from mpl_toolkits.axes_grid1 import make_axes_locatable
+import argparse
# Parameters users should/may tweak
# plot all groups and all photon quantities
plot_all_data = True
-# snapshot basename
-snapshot_base = "output"
# fancy up the plots a bit?
fancy = True
@@ -49,15 +48,20 @@ projection_kwargs = {"resolution": 512, "parallel": True}
# -----------------------------------------------------------------------
-
+mpl.rcParams["text.usetex"] = True
# Read in cmdline arg: Are we plotting only one snapshot, or all?
plot_all = False
-try:
- snapnr = int(sys.argv[1])
-except IndexError:
- plot_all = True
-mpl.rcParams["text.usetex"] = True
+
+def parse_args():
+ parser = argparse.ArgumentParser()
+ parser.add_argument(
+ "-n", "--snapshot-number", help="Number of snapshot to plot", type=int
+ )
+ parser.add_argument("-b", "--basename", help="Snapshot basename", default="output")
+
+ args = parser.parse_args()
+ return args
def get_snapshot_list(snapshot_basename="output"):
@@ -75,6 +79,9 @@ def get_snapshot_list(snapshot_basename="output"):
snaplist.append(f)
snaplist = sorted(snaplist)
+ if len(snaplist) == 0:
+ print(f"No snapshots with base '{snapshot_basename}' found")
+ sys.exit(1)
else:
fname = snapshot_basename + "_" + str(snapnr).zfill(4) + ".hdf5"
@@ -224,7 +231,16 @@ def plot_ionization(filename):
if __name__ == "__main__":
+ # Get command line arguments
+ args = parse_args()
+
+ if args.snapshot_number:
+ plot_all = False
+ snapnr = int(args.snapshot_numbeR)
+ else:
+ plot_all = True
+ snapshot_base = args.basename
snaplist = get_snapshot_list(snapshot_base)
for f in snaplist:
diff --git a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/runCosmo.sh b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/runCosmo.sh
new file mode 100755
index 0000000000000000000000000000000000000000..b3e9edd3c15c01d01336fb05e912201de32a1678
--- /dev/null
+++ b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/runCosmo.sh
@@ -0,0 +1,30 @@
+#!/bin/bash
+
+# make run.sh fail if a subcommand fails
+set -e
+set -o pipefail
+
+if [ ! -e glassPlane_64.hdf5 ]
+then
+ echo "Fetching initial glass file ..."
+ ./getGlass.sh
+fi
+
+if [ ! -f 'advect_ions.hdf5' ]; then
+ echo "Generating ICs"
+ python3 makeIC.py
+fi
+
+# Run SWIFT with RT
+../../../swift \
+ --hydro \
+ --cosmology \
+ --threads=4 \
+ --stars \
+ --external-gravity \
+ --feedback \
+ --radiation \
+ advect_ions_cosmo.yml 2>&1 | tee output.log
+
+python3 plotIonization.py -b cosmo
+python3 testIonization.py -b cosmo
diff --git a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/snapshot_times_high_redshift b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/snapshot_times_high_redshift
new file mode 100644
index 0000000000000000000000000000000000000000..e504014177837a7b2c1d134d1c4a0d5a785d0044
--- /dev/null
+++ b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/snapshot_times_high_redshift
@@ -0,0 +1,12 @@
+# Redshift
+99.98460572382712
+96.47031888378449
+93.23367859034735
+90.25299848922855
+87.49735694548497
+84.9407463291136
+82.56112036003019
+80.33965676782498
+78.26018022924342
+76.30870609169932
+74.47307621234665
diff --git a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/testIonization.py b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/testIonization.py
index 3c737462e2ea28f61897c34c6b0d7642626aa463..da002bea6ecbb36e819af6a71e12cc452a0dfa48 100755
--- a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/testIonization.py
+++ b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/testIonization.py
@@ -25,11 +25,19 @@
# ----------------------------------------------------
import os
+import sys
import swiftsimio
from matplotlib import pyplot as plt
+import argparse
-snapshot_base = "output"
+
+def parse_args():
+ parser = argparse.ArgumentParser()
+ parser.add_argument("-b", "--basename", help="Snapshot basename", default="output")
+
+ args = parser.parse_args()
+ return args
def get_snapshot_list(snapshot_basename="output"):
@@ -46,6 +54,9 @@ def get_snapshot_list(snapshot_basename="output"):
snaplist.append(f)
snaplist = sorted(snaplist)
+ if len(snaplist) == 0:
+ print(f"No snapshots with base '{snapshot_basename}' found")
+ sys.exit(1)
return snaplist
@@ -61,6 +72,11 @@ def compare_data(snaplist):
HeII = []
HeIII = []
+ if "cosmo" in snaplist[0]:
+ savename = "total_abundancies_cosmo.png"
+ else:
+ savename = "total_abundancies.png"
+
for filename in snaplist:
data = swiftsimio.load(filename)
@@ -86,12 +102,15 @@ def compare_data(snaplist):
# plt.show()
plt.tight_layout()
- plt.savefig("total_abundancies.png", dpi=200)
+ plt.savefig(savename, dpi=200)
return
if __name__ == "__main__":
+ # Get command line arguments
+ args = parse_args()
+ snapshot_base = args.basename
snaplist = get_snapshot_list(snapshot_base)
compare_data(snaplist)