Skip to content
Snippets Groups Projects
Commit 00e9de0a authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'GEARRT_cosmo_IonMassFraction_example' into 'master'

Cosmology IonMassFraction example GEAR-RT

See merge request !1857
parents d964bb05 60deb638
No related branches found
No related tags found
2 merge requests!1891Merge master into Zoom merge,!1857Cosmology IonMassFraction example GEAR-RT
...@@ -13,7 +13,7 @@ InternalUnitSystem: ...@@ -13,7 +13,7 @@ InternalUnitSystem:
TimeIntegration: TimeIntegration:
max_nr_rt_subcycles: 1 max_nr_rt_subcycles: 1
time_begin: 0. # The starting time of the simulation (in internal units). 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_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). dt_max: 1.e-03 # The maximal time-step size of the simulation (in internal units).
...@@ -22,7 +22,7 @@ TimeIntegration: ...@@ -22,7 +22,7 @@ TimeIntegration:
Snapshots: Snapshots:
basename: output # Common part of the name of output files basename: output # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units) 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 # Parameters governing the conserved quantities statistics
Statistics: Statistics:
......
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
...@@ -43,7 +43,7 @@ if __name__ == "__main__": ...@@ -43,7 +43,7 @@ if __name__ == "__main__":
# Set up metadata # Set up metadata
unitL = unyt.Mpc unitL = unyt.Mpc
edgelen = 2 * 15 * 1e-3 * unitL # 30 kpc edgelen = 1 * unyt.Mpc
edgelen = edgelen.to(unitL) edgelen = edgelen.to(unitL)
boxsize = np.array([1.0, 1.0, 0.0]) * edgelen boxsize = np.array([1.0, 1.0, 0.0]) * edgelen
......
...@@ -31,13 +31,12 @@ import swiftsimio ...@@ -31,13 +31,12 @@ import swiftsimio
from matplotlib import pyplot as plt from matplotlib import pyplot as plt
from matplotlib.colors import SymLogNorm from matplotlib.colors import SymLogNorm
from mpl_toolkits.axes_grid1 import make_axes_locatable from mpl_toolkits.axes_grid1 import make_axes_locatable
import argparse
# Parameters users should/may tweak # Parameters users should/may tweak
# plot all groups and all photon quantities # plot all groups and all photon quantities
plot_all_data = True plot_all_data = True
# snapshot basename
snapshot_base = "output"
# fancy up the plots a bit? # fancy up the plots a bit?
fancy = True fancy = True
...@@ -49,15 +48,20 @@ projection_kwargs = {"resolution": 512, "parallel": 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? # Read in cmdline arg: Are we plotting only one snapshot, or all?
plot_all = False 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"): def get_snapshot_list(snapshot_basename="output"):
...@@ -75,6 +79,9 @@ def get_snapshot_list(snapshot_basename="output"): ...@@ -75,6 +79,9 @@ def get_snapshot_list(snapshot_basename="output"):
snaplist.append(f) snaplist.append(f)
snaplist = sorted(snaplist) snaplist = sorted(snaplist)
if len(snaplist) == 0:
print(f"No snapshots with base '{snapshot_basename}' found")
sys.exit(1)
else: else:
fname = snapshot_basename + "_" + str(snapnr).zfill(4) + ".hdf5" fname = snapshot_basename + "_" + str(snapnr).zfill(4) + ".hdf5"
...@@ -224,7 +231,16 @@ def plot_ionization(filename): ...@@ -224,7 +231,16 @@ def plot_ionization(filename):
if __name__ == "__main__": 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) snaplist = get_snapshot_list(snapshot_base)
for f in snaplist: for f in snaplist:
......
#!/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
# Redshift
99.98460572382712
96.47031888378449
93.23367859034735
90.25299848922855
87.49735694548497
84.9407463291136
82.56112036003019
80.33965676782498
78.26018022924342
76.30870609169932
74.47307621234665
...@@ -25,11 +25,19 @@ ...@@ -25,11 +25,19 @@
# ---------------------------------------------------- # ----------------------------------------------------
import os import os
import sys
import swiftsimio import swiftsimio
from matplotlib import pyplot as plt 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"): def get_snapshot_list(snapshot_basename="output"):
...@@ -46,6 +54,9 @@ def get_snapshot_list(snapshot_basename="output"): ...@@ -46,6 +54,9 @@ def get_snapshot_list(snapshot_basename="output"):
snaplist.append(f) snaplist.append(f)
snaplist = sorted(snaplist) snaplist = sorted(snaplist)
if len(snaplist) == 0:
print(f"No snapshots with base '{snapshot_basename}' found")
sys.exit(1)
return snaplist return snaplist
...@@ -61,6 +72,11 @@ def compare_data(snaplist): ...@@ -61,6 +72,11 @@ def compare_data(snaplist):
HeII = [] HeII = []
HeIII = [] HeIII = []
if "cosmo" in snaplist[0]:
savename = "total_abundancies_cosmo.png"
else:
savename = "total_abundancies.png"
for filename in snaplist: for filename in snaplist:
data = swiftsimio.load(filename) data = swiftsimio.load(filename)
...@@ -86,12 +102,15 @@ def compare_data(snaplist): ...@@ -86,12 +102,15 @@ def compare_data(snaplist):
# plt.show() # plt.show()
plt.tight_layout() plt.tight_layout()
plt.savefig("total_abundancies.png", dpi=200) plt.savefig(savename, dpi=200)
return return
if __name__ == "__main__": if __name__ == "__main__":
# Get command line arguments
args = parse_args()
snapshot_base = args.basename
snaplist = get_snapshot_list(snapshot_base) snaplist = get_snapshot_list(snapshot_base)
compare_data(snaplist) compare_data(snaplist)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment