diff --git a/examples/RadiativeTransferTests/Advection_1D/makeIC.py b/examples/RadiativeTransferTests/Advection_1D/makeIC.py
index 508b96647d8360a76368090d14186e6e0c07d808..4b2abcb4a1563f59179eec151327e6135bb4cb7c 100755
--- a/examples/RadiativeTransferTests/Advection_1D/makeIC.py
+++ b/examples/RadiativeTransferTests/Advection_1D/makeIC.py
@@ -35,7 +35,6 @@ from swiftsimio import Writer
import unyt
import numpy as np
import h5py
-from matplotlib import pyplot as plt
# define unit system to use
unitsystem = unyt.unit_systems.cgs_unit_system
@@ -83,7 +82,7 @@ def initial_condition(x):
# Assuming all photons flow in only one direction
# (optically thin regime, "free streaming limit"),
# we have that |F| = c * E
- F = np.zeros(3, dtype=np.float32)
+ F = np.zeros(3, dtype=np.float64)
F[0] = unyt.c.to(unitsystem["length"] / unitsystem["time"]) * E
E_list.append(E)
@@ -99,7 +98,7 @@ def initial_condition(x):
else:
E = 1.0
- F = np.zeros(3, dtype=np.float32)
+ F = np.zeros(3, dtype=np.float64)
F[0] = unyt.c.to(unitsystem["length"] / unitsystem["time"]) * E
E_list.append(E)
@@ -112,7 +111,7 @@ def initial_condition(x):
amplitude = 2.0
E = amplitude * np.exp(-(x[0] - mean) ** 2 / (2 * sigma ** 2))
- F = np.zeros(3, dtype=np.float32)
+ F = np.zeros(3, dtype=np.float64)
F[0] = unyt.c.to(unitsystem["length"] / unitsystem["time"]) * E
E_list.append(E)
@@ -123,7 +122,7 @@ def initial_condition(x):
if __name__ == "__main__":
- xp = unyt.unyt_array(np.zeros((n_p, 3), dtype=np.float32), boxsize.units)
+ xp = unyt.unyt_array(np.zeros((n_p, 3), dtype=np.float64), boxsize.units)
dx = boxsize / n_p
@@ -134,9 +133,9 @@ if __name__ == "__main__":
w.gas.coordinates = xp
w.gas.velocities = np.zeros(xp.shape) * (unyt.cm / unyt.s)
- w.gas.masses = np.ones(xp.shape[0], dtype=np.float32) * 1000 * unyt.g
+ w.gas.masses = np.ones(xp.shape[0], dtype=np.float64) * 1000 * unyt.g
w.gas.internal_energy = (
- np.ones(xp.shape[0], dtype=np.float32) * (300.0 * unyt.kb * unyt.K) / (unyt.g)
+ np.ones(xp.shape[0], dtype=np.float64) * (300.0 * unyt.kb * unyt.K) / (unyt.g)
)
# Generate initial guess for smoothing lengths based on MIPS
@@ -172,6 +171,7 @@ if __name__ == "__main__":
Fsetname = "PhotonFluxesGroup{0:d}".format(g + 1)
parts[Fsetname][p] = Flux[g]
+ # from matplotlib import pyplot as plt
# plt.figure()
# for g in range(nPhotonGroups):
# # Esetname = "PhotonEnergiesGroup{0:d}".format(g+1)
diff --git a/examples/RadiativeTransferTests/Advection_1D/plotSolution.py b/examples/RadiativeTransferTests/Advection_1D/plotSolution.py
index 00becd2d8e2ac96df2b9299c90fb9ce0b3ba6d82..8209e6fc58226a84741fe2854282c063e570545b 100755
--- a/examples/RadiativeTransferTests/Advection_1D/plotSolution.py
+++ b/examples/RadiativeTransferTests/Advection_1D/plotSolution.py
@@ -148,7 +148,7 @@ def plot_photons(filename, energy_boundaries=None, flux_boundaries=None):
advected_positions[overshooters] -= boxsize
analytical_solutions = np.zeros(
- (part_positions.shape[0], ngroups), dtype=np.float
+ (part_positions.shape[0], ngroups), dtype=np.float64
)
for p in range(part_positions.shape[0]):
E, F = initial_condition(advected_positions[p])
@@ -319,7 +319,6 @@ def get_minmax_vals(snaplist):
for direction in ["X"]:
# for direction in ["X", "Y", "Z"]:
- new_attribute_str = "radiation_flux" + str(g + 1) + direction
f = getattr(data.gas.photon_fluxes, "Group" + str(g + 1) + direction)
fluxmin_group.append(f.min())
fluxmax_group.append(f.max())
diff --git a/examples/RadiativeTransferTests/Advection_2D/makeIC.py b/examples/RadiativeTransferTests/Advection_2D/makeIC.py
index 5e2d3507ac39d3176cdae469c4de7ad50b734e6a..4273e2d9da423b9b0b5c4e883f48a0c52bef3655 100755
--- a/examples/RadiativeTransferTests/Advection_2D/makeIC.py
+++ b/examples/RadiativeTransferTests/Advection_2D/makeIC.py
@@ -35,7 +35,6 @@ from swiftsimio import Writer
import unyt
import numpy as np
import h5py
-from matplotlib import pyplot as plt
# define unit system to use
@@ -82,7 +81,7 @@ def initial_condition(x):
# Assuming all photons flow in only one direction
# (optically thin regime, "free streaming limit"),
# we have that |F| = c * E
- F = np.zeros(3, dtype=np.float32)
+ F = np.zeros(3, dtype=np.float64)
F[0] = c * E
E_list.append(E)
@@ -98,7 +97,7 @@ def initial_condition(x):
else:
E = 1.0
- F = np.zeros(3, dtype=np.float32)
+ F = np.zeros(3, dtype=np.float64)
F[1] = c * E
E_list.append(E)
@@ -116,7 +115,7 @@ def initial_condition(x):
* np.exp(-((x[0] - mean) ** 2 + (x[1] - mean) ** 2) / (2 * sigma ** 2))
+ baseline
)
- F = np.zeros(3, dtype=np.float32)
+ F = np.zeros(3, dtype=np.float64)
F[0] = c * E * 1.414213562 # sqrt(2)
F[1] = c * E * 1.414213562 # sqrt(2)
@@ -135,13 +134,13 @@ def initial_condition(x):
unit_vector = (dx / r, dy / r)
E = 1.0
- F = np.zeros(3, dtype=np.float32)
+ F = np.zeros(3, dtype=np.float64)
F[0] = unit_vector[0] * c * E
F[1] = unit_vector[1] * c * E
else:
E = 0.0
- F = np.zeros(3, dtype=np.float32)
+ F = np.zeros(3, dtype=np.float64)
E_list.append(E)
F_list.append(F)
@@ -166,9 +165,9 @@ if __name__ == "__main__":
w.gas.coordinates = pos
w.gas.velocities = np.zeros((numPart, 3)) * (unyt.cm / unyt.s)
- w.gas.masses = np.ones(numPart, dtype=np.float32) * 1000 * unyt.g
+ w.gas.masses = np.ones(numPart, dtype=np.float64) * 1000 * unyt.g
w.gas.internal_energy = (
- np.ones(numPart, dtype=np.float32) * (300.0 * unyt.kb * unyt.K) / (unyt.g)
+ np.ones(numPart, dtype=np.float64) * (300.0 * unyt.kb * unyt.K) / (unyt.g)
)
# Generate initial guess for smoothing lengths based on MIPS
diff --git a/examples/RadiativeTransferTests/Advection_2D/plotSolution.py b/examples/RadiativeTransferTests/Advection_2D/plotSolution.py
index 433aa76160716a77a6e7e21c653df1e8fad4d62c..29bb0b012fe75c6eaddcd1cb0ffcbed92e6c2f39 100755
--- a/examples/RadiativeTransferTests/Advection_2D/plotSolution.py
+++ b/examples/RadiativeTransferTests/Advection_2D/plotSolution.py
@@ -29,7 +29,6 @@
import sys
import os
import swiftsimio
-import numpy as np
import gc
from matplotlib import pyplot as plt
import matplotlib as mpl
@@ -133,7 +132,7 @@ def plot_photons(filename, energy_boundaries=None, flux_boundaries=None):
f *= data.gas.masses
setattr(data.gas, new_attribute_str, f)
- # get mass surface density projection that we'll use to remove density dependence in impage
+ # get mass surface density projection that we'll use to remove density dependence in image
mass_map = swiftsimio.visualisation.projection.project_gas(
data, project="masses", **projection_kwargs
)
@@ -271,7 +270,6 @@ def get_minmax_vals(snaplist):
dirmin = []
dirmax = []
for direction in ["X", "Y"]:
- new_attribute_str = "radiation_flux" + str(g + 1) + direction
f = getattr(data.gas.photon_fluxes, "Group" + str(g + 1) + direction)
dirmin.append(f.min())
dirmax.append(f.max())
diff --git a/examples/RadiativeTransferTests/Advection_2D/plotSolutionScatter.py b/examples/RadiativeTransferTests/Advection_2D/plotSolutionScatter.py
index 3e8491ec5d709f11c271d7a9da6d0eb87f3270f9..a5180c5507a90573b1f6e60a40782b2b523af81b 100755
--- a/examples/RadiativeTransferTests/Advection_2D/plotSolutionScatter.py
+++ b/examples/RadiativeTransferTests/Advection_2D/plotSolutionScatter.py
@@ -29,7 +29,6 @@
import sys
import os
import swiftsimio
-import numpy as np
import gc
from matplotlib import pyplot as plt
import matplotlib as mpl
@@ -87,22 +86,11 @@ def get_snapshot_list(snapshot_basename="output"):
return snaplist
-def set_colorbar(ax, im):
- divider = make_axes_locatable(ax)
- cax = divider.append_axes("right", size="5%", pad=0.05)
- plt.colorbar(im, cax=cax)
- return
-
-
-def plot_photons(filename, energy_boundaries=None, flux_boundaries=None):
+def plot_photons(filename):
"""
Create the actual plot.
filename: file to work with
- energy_boundaries: list of [E_min, E_max] for each photon group.
- If none, limits are set automatically.
- flux_boundaries: list of [F_min, F_max] for each photon group.
- If none, limits are set automatically.
"""
print("working on", filename)
@@ -124,7 +112,6 @@ def plot_photons(filename, energy_boundaries=None, flux_boundaries=None):
if plot_all_data:
# prepare also the fluxes
- # for direction in ["X", "Y", "Z"]:
for direction in ["X", "Y"]:
new_attribute_str = "radiation_flux" + str(g + 1) + direction
f = getattr(data.gas.photon_fluxes, "Group" + str(g + 1) + direction)
@@ -204,78 +191,9 @@ def plot_photons(filename, energy_boundaries=None, flux_boundaries=None):
return
-def get_minmax_vals(snaplist):
- """
- Find minimal and maximal values for energy and flux
- so you can fix axes limits over all snapshots
-
- snaplist: list of snapshot filenames
-
- returns:
-
- energy_boundaries: list of [E_min, E_max] for each photon group
- flux_boundaries: list of [Fx_min, Fy_max] for each photon group
- """
-
- emins = []
- emaxs = []
- fmins = []
- fmaxs = []
-
- for filename in snaplist:
-
- data = swiftsimio.load(filename)
- meta = data.metadata
-
- ngroups = int(meta.subgrid_scheme["PhotonGroupNumber"])
- emin_group = []
- emax_group = []
- fluxmin_group = []
- fluxmax_group = []
-
- for g in range(ngroups):
- en = getattr(data.gas.photon_energies, "group" + str(g + 1))
- emin_group.append(en.min())
- emax_group.append(en.max())
-
- dirmin = []
- dirmax = []
- for direction in ["X", "Y"]:
- new_attribute_str = "radiation_flux" + str(g + 1) + direction
- f = getattr(data.gas.photon_fluxes, "Group" + str(g + 1) + direction)
- dirmin.append(f.min())
- dirmax.append(f.max())
- fluxmin_group.append(min(dirmin))
- fluxmax_group.append(max(dirmax))
-
- emins.append(emin_group)
- emaxs.append(emax_group)
- fmins.append(fluxmin_group)
- fmaxs.append(fluxmax_group)
-
- energy_boundaries = []
- flux_boundaries = []
- for g in range(ngroups):
- emin = min([emins[f][g] for f in range(len(snaplist))])
- emax = max([emaxs[f][g] for f in range(len(snaplist))])
- energy_boundaries.append([emin, emax])
- fmin = min([fmins[f][g] for f in range(len(snaplist))])
- fmax = max([fmaxs[f][g] for f in range(len(snaplist))])
- flux_boundaries.append([fmin, fmax])
-
- return energy_boundaries, flux_boundaries
-
-
if __name__ == "__main__":
snaplist = get_snapshot_list(snapshot_base)
- if fancy:
- energy_boundaries, flux_boundaries = get_minmax_vals(snaplist)
- else:
- energy_boundaries = None
- flux_boundaries = None
for f in snaplist:
- plot_photons(
- f, energy_boundaries=energy_boundaries, flux_boundaries=flux_boundaries
- )
+ plot_photons(f)
diff --git a/examples/RadiativeTransferTests/Advection_2D/rt_advection2D.yml b/examples/RadiativeTransferTests/Advection_2D/rt_advection2D.yml
index 8c1b9d88f980109b1446960d0741f13015562e82..5838bc59d8a8405ec7b6f7eefd36ba1c8ddf4635 100644
--- a/examples/RadiativeTransferTests/Advection_2D/rt_advection2D.yml
+++ b/examples/RadiativeTransferTests/Advection_2D/rt_advection2D.yml
@@ -29,7 +29,7 @@ Statistics:
# Parameters for the hydrodynamics scheme
SPH:
- resolution_eta: 3.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ 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.6 # Courant-Friedrich-Levy condition for time integration.
minimal_temperature: 10. # Kelvin
diff --git a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/README b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/README
new file mode 100644
index 0000000000000000000000000000000000000000..0f0df76c4d5ffd7e4b8bc79fe5f2e2d4c8e7ce11
--- /dev/null
+++ b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/README
@@ -0,0 +1,9 @@
+Ion Mass Fraction Advection Test
+
+In some schemes, like the meshless (gizmo) scheme, particles exchange
+masses via fluxes. This example tests whether the mass flux exchange
+is done correctly by setting up regions of special ionization states
+and advecting them through time.
+
+To use the GEAR RT model, compile with :
+ --with-rt=GEAR_1 --with-rt-riemann-solver=GLF --with-hydro-dimension=2 --with-hydro=gizmo-mfv --with-riemann-solver=hllc --with-stars=GEAR --with-feedback=none --with-grackle=$GRACKLE_ROOT
diff --git a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/advect_ions.yml b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/advect_ions.yml
new file mode 100644
index 0000000000000000000000000000000000000000..9db69ab1cefce55353de7b72580a5e692cf707b3
--- /dev/null
+++ b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/advect_ions.yml
@@ -0,0 +1,57 @@
+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:
+ 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: 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
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ time_first: 0.
+ delta_time: 5e-4 # 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.
+
+Scheduler:
+ max_top_level_cells: 12
+
+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
+ photon_groups_Hz: [3.288e15, 5.945e15, 13.157e15] # Photon frequency group bin edges in Hz. Needs to be 1 less than the number of groups (N) requested during the configuration (--with-RT=GEAR_N). Outer edges of zero and infinity are assumed.
+ use_const_emission_rates: 1 # (Optional) use constant emission rates for stars as defined with star_emission_rates_erg_LSol parameter
+ star_emission_rates_LSol: [0., 0., 0., 0.] # (Optional) constant star emission rates for each photon frequency group to use if use_constant_emission_rates is set, in units of Solar Luminosity.
+ hydrogen_mass_fraction: 0.50 # 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
+
diff --git a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/getGlass.sh b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..a958311b6de07663c7f7c70cff9e95d86ce3efb2
--- /dev/null
+++ b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/getGlass.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassPlane_64.hdf5
diff --git a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/makeIC.py b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/makeIC.py
new file mode 100755
index 0000000000000000000000000000000000000000..f9f0c1afcd6ac9092e0d3c7f50d0bfa13198a0ae
--- /dev/null
+++ b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/makeIC.py
@@ -0,0 +1,97 @@
+#!/usr/bin/env python3
+
+# ---------------------------------------------------------------------
+# Test the diffusion/advection of ions by creating regions with/without
+# initial ions. Run without radiation.
+# ---------------------------------------------------------------------
+
+from swiftsimio import Writer
+from swiftsimio.units import cosmo_units
+import unyt
+import numpy as np
+import h5py
+
+gamma = 5.0 / 3.0
+outputfilename = "advect_ions.hdf5"
+
+if __name__ == "__main__":
+
+ glass = h5py.File("glassPlane_64.hdf5", "r")
+ parts = glass["PartType0"]
+ xp = parts["Coordinates"][:]
+ h = parts["SmoothingLength"][:]
+ glass.close()
+
+ # Set up metadata
+ unitL = unyt.Mpc
+ edgelen = 2 * 15 * 1e-3 * unitL # 30 kpc
+ edgelen = edgelen.to(unitL)
+ boxsize = np.array([1.0, 1.0, 0.0]) * edgelen
+
+ xp *= edgelen
+ h *= edgelen
+
+ w = Writer(unit_system=cosmo_units, box_size=boxsize, dimension=2)
+
+ # write particle positions ans smoothing lengths
+ w.gas.coordinates = xp
+ w.gas.smoothing_length = h
+
+ # get gas masses
+ mpart = 1.6e5 * unyt.Msun
+ masses = np.ones(xp.shape[0], dtype=np.float64) * mpart
+ # change some gas masses
+ mask = xp[:, 0] > 0.5 * edgelen
+ masses[mask] *= 3
+ w.gas.masses = masses
+
+ w.gas.internal_energy = (
+ np.ones(xp.shape[0], dtype=np.float64) * 1.25e6 * unyt.m ** 2 / unyt.s ** 2
+ )
+
+ # get velocities
+ vels = np.zeros((xp.shape[0], 3))
+ vels[:, 0] = -1.0
+ vels[:, 1] = +1.0
+ w.gas.velocities = vels * 1000 * cosmo_units["length"] / cosmo_units["time"]
+
+ w.write(outputfilename)
+
+ # Now open file back up again and add RT data.
+ F = h5py.File(outputfilename, "r+")
+ header = F["Header"]
+ nparts = header.attrs["NumPart_ThisFile"][0]
+ parts = F["/PartType0"]
+ pos = parts["Coordinates"]
+
+ # Create initial ionization species mass fractions.
+ HIdata = np.ones((nparts), dtype=np.float32) * 0.76
+ HIIdata = np.zeros((nparts), dtype=np.float32)
+ HeIdata = np.ones((nparts), dtype=np.float32) * 0.24
+ HeIIdata = np.zeros((nparts), dtype=np.float32)
+ HeIIIdata = np.zeros((nparts), dtype=np.float32)
+
+ mask1 = pos[:, 0] > edgelen / 3
+ mask1 = np.logical_and(mask1, pos[:, 0] < 2 * edgelen / 3)
+ HIdata[mask1] = 0.26
+ HIIdata[mask1] = 0.5
+
+ mask2 = pos[:, 1] > edgelen / 4
+ mask2 = np.logical_and(mask2, pos[:, 1] < edgelen / 2)
+ HeIdata[mask2] = 0.1
+ HeIIdata[mask2] = 0.14
+
+ mask3 = pos[:, 1] > edgelen / 2
+ mask3 = np.logical_and(mask3, pos[:, 1] < 3 * edgelen / 4)
+ HeIdata[mask3] = 0.05
+ HeIIdata[mask3] = 0.09
+ HeIIIdata[mask3] = 0.1
+
+ parts.create_dataset("MassFractionHI", data=HIdata)
+ parts.create_dataset("MassFractionHII", data=HIIdata)
+ parts.create_dataset("MassFractionHeI", data=HeIdata)
+ parts.create_dataset("MassFractionHeII", data=HeIIdata)
+ parts.create_dataset("MassFractionHeIII", data=HeIIIdata)
+
+ # close up, and we're done!
+ F.close()
diff --git a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/plotIonization.py b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/plotIonization.py
new file mode 100755
index 0000000000000000000000000000000000000000..6f0112bd09d1f33ccf145da888c112a2ccdd45d1
--- /dev/null
+++ b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/plotIonization.py
@@ -0,0 +1,214 @@
+#!/usr/bin/env python3
+
+# ----------------------------------------------------
+# Plot the ionizing species
+# ----------------------------------------------------
+
+import sys
+import os
+import swiftsimio
+import numpy as np
+import gc
+import unyt
+from matplotlib import pyplot as plt
+import matplotlib as mpl
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+from matplotlib.colors import SymLogNorm
+
+# 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
+
+# parameters for imshow plots
+imshow_kwargs = {"origin": "lower", "cmap": "brg"}
+
+# parameters for swiftsimio projections
+projection_kwargs = {"resolution": 512, "parallel": 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 get_snapshot_list(snapshot_basename="output"):
+ """
+ Find the snapshot(s) that are to be plotted
+ and return their names as list
+ """
+
+ snaplist = []
+
+ if plot_all:
+ dirlist = os.listdir()
+ for f in dirlist:
+ if f.startswith(snapshot_basename) and f.endswith("hdf5"):
+ snaplist.append(f)
+
+ snaplist = sorted(snaplist)
+
+ else:
+ fname = snapshot_basename + "_" + str(snapnr).zfill(4) + ".hdf5"
+ if not os.path.exists(fname):
+ print("Didn't find file", fname)
+ quit(1)
+ snaplist.append(fname)
+
+ return snaplist
+
+
+def set_colorbar(ax, im):
+ """
+ Adapt the colorbar a bit for axis object <ax> and
+ imshow instance <im>
+ """
+ divider = make_axes_locatable(ax)
+ cax = divider.append_axes("right", size="5%", pad=0.05)
+ plt.colorbar(im, cax=cax)
+ return
+
+
+def plot_ionization(filename):
+ """
+ Create and save the plot
+ """
+ print("working on", filename)
+
+ data = swiftsimio.load(filename)
+ meta = data.metadata
+
+ global imshow_kwargs
+ imshow_kwargs["extent"] = [
+ 0.01 * meta.boxsize[0].v,
+ 0.99 * meta.boxsize[0].v,
+ 0.01 * meta.boxsize[1].v,
+ 0.99 * meta.boxsize[1].v,
+ ]
+
+ mass_map = swiftsimio.visualisation.projection.project_gas(
+ data, project="masses", **projection_kwargs
+ )
+
+ data.gas.mXHI = data.gas.ion_mass_fractions.HI * data.gas.masses
+ data.gas.mXHII = data.gas.ion_mass_fractions.HII * data.gas.masses
+ data.gas.mXHeI = data.gas.ion_mass_fractions.HeI * data.gas.masses
+ data.gas.mXHeII = data.gas.ion_mass_fractions.HeII * data.gas.masses
+ data.gas.mXHeIII = data.gas.ion_mass_fractions.HeIII * data.gas.masses
+
+ mass_weighted_XHImap = swiftsimio.visualisation.projection.project_gas(
+ data, project="mXHI", **projection_kwargs
+ )
+ mass_weighted_XHIImap = swiftsimio.visualisation.projection.project_gas(
+ data, project="mXHII", **projection_kwargs
+ )
+ mass_weighted_XHeImap = swiftsimio.visualisation.projection.project_gas(
+ data, project="mXHeI", **projection_kwargs
+ )
+ mass_weighted_XHeIImap = swiftsimio.visualisation.projection.project_gas(
+ data, project="mXHeII", **projection_kwargs
+ )
+ mass_weighted_XHeIIImap = swiftsimio.visualisation.projection.project_gas(
+ data, project="mXHeIII", **projection_kwargs
+ )
+
+ XHImap = mass_weighted_XHImap / mass_map
+ XHIImap = mass_weighted_XHIImap / mass_map
+ XHeImap = mass_weighted_XHeImap / mass_map
+ XHeIImap = mass_weighted_XHeIImap / mass_map
+ XHeIIImap = mass_weighted_XHeIIImap / mass_map
+
+ fig = plt.figure(figsize=(20, 8), dpi=200)
+ figname = filename[:-5] + "-mass-fractions.png"
+
+ ax1 = fig.add_subplot(251)
+ ax2 = fig.add_subplot(252)
+ ax3 = fig.add_subplot(253)
+ ax4 = fig.add_subplot(254)
+ ax5 = fig.add_subplot(255)
+ ax6 = fig.add_subplot(256)
+ ax7 = fig.add_subplot(257)
+ ax8 = fig.add_subplot(258)
+ ax9 = fig.add_subplot(259)
+ ax10 = fig.add_subplot(2, 5, 10)
+
+ # im0 = ax0.imshow(mass_map.T, **imshow_kwargs)
+ # set_colorbar(ax0, im0)
+ # ax0.set_title("Mass")
+
+ imshow_kwargs_ions = imshow_kwargs.copy()
+ imshow_kwargs_ions["norm"] = SymLogNorm(vmin=0.0, vmax=1.0, linthresh=1e-3, base=10)
+
+ im1 = ax1.imshow(XHImap.T, **imshow_kwargs_ions)
+ set_colorbar(ax1, im1)
+ ax1.set_title("HI Mass Fraction")
+
+ im2 = ax2.imshow(XHIImap.T, **imshow_kwargs_ions)
+ set_colorbar(ax2, im2)
+ ax2.set_title("HII Mass Fraction")
+
+ im3 = ax3.imshow(XHeImap.T, **imshow_kwargs_ions)
+ set_colorbar(ax3, im3)
+ ax3.set_title("HeI Mass Fraction")
+
+ im4 = ax4.imshow(XHeIImap.T, **imshow_kwargs_ions)
+ set_colorbar(ax4, im4)
+ ax4.set_title("HeII Mass Fraction")
+
+ im5 = ax5.imshow(XHeIIImap.T, **imshow_kwargs_ions)
+ set_colorbar(ax5, im5)
+ ax5.set_title("HeIII Mass Fraction")
+
+ im6 = ax6.imshow(mass_weighted_XHImap.T, **imshow_kwargs)
+ set_colorbar(ax6, im6)
+ ax6.set_title("HI Mass")
+
+ im7 = ax7.imshow(mass_weighted_XHIImap.T, **imshow_kwargs)
+ set_colorbar(ax7, im7)
+ ax7.set_title("HII Mass")
+
+ im8 = ax8.imshow(mass_weighted_XHeImap.T, **imshow_kwargs)
+ set_colorbar(ax8, im8)
+ ax8.set_title("HeI Mass")
+
+ im9 = ax9.imshow(mass_weighted_XHeIImap.T, **imshow_kwargs)
+ set_colorbar(ax9, im9)
+ ax9.set_title("HeII Mass")
+
+ im10 = ax10.imshow(mass_weighted_XHeIIImap.T, **imshow_kwargs)
+ set_colorbar(ax10, im10)
+ ax10.set_title("HeIII Mass")
+
+ for ax in [ax1, ax2, ax3, ax4, ax5, ax6, ax7, ax8, ax9, ax10]:
+ ax.set_xlabel("[kpc]")
+ ax.set_ylabel("[kpc]")
+
+ title = filename.replace("_", "\_") # exception handle underscore for latex
+ if meta.cosmology is not None:
+ title += ", $z$ = {0:.2e}".format(meta.z)
+ title += ", $t$ = {0:.2e}".format(meta.time.to("Myr"))
+ fig.suptitle(title)
+
+ plt.tight_layout()
+ plt.savefig(figname)
+ plt.close()
+ return
+
+
+if __name__ == "__main__":
+
+ snaplist = get_snapshot_list(snapshot_base)
+
+ for f in snaplist:
+ plot_ionization(f)
diff --git a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/run.sh b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..51346e2509df292d3368b7bcbd08b982d626ad0d
--- /dev/null
+++ b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/run.sh
@@ -0,0 +1,25 @@
+#!/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 --threads=4 --stars --external-gravity \
+ --feedback --radiation \
+ advect_ions.yml 2>&1 | tee output.log
+
+python3 plotIonization.py
+python3 testIonization.py
diff --git a/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/testIonization.py b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/testIonization.py
new file mode 100755
index 0000000000000000000000000000000000000000..fa2cd1fea1a6414f095cf4f800f1d339c86ae762
--- /dev/null
+++ b/examples/RadiativeTransferTests/IonMassFractionAdvectionTest_2D/testIonization.py
@@ -0,0 +1,84 @@
+#!/usr/bin/env python3
+
+# ----------------------------------------------------
+# Check that the total amount of ionized species
+# remains constant
+# ----------------------------------------------------
+
+import sys
+import os
+import swiftsimio
+import numpy as np
+import gc
+import unyt
+from matplotlib import pyplot as plt
+import matplotlib as mpl
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+from matplotlib.colors import SymLogNorm
+
+snapshot_base = "output"
+
+
+def get_snapshot_list(snapshot_basename="output"):
+ """
+ Find the snapshot(s) that are to be plotted
+ and return their names as list
+ """
+
+ snaplist = []
+
+ dirlist = os.listdir()
+ for f in dirlist:
+ if f.startswith(snapshot_basename) and f.endswith("hdf5"):
+ snaplist.append(f)
+
+ snaplist = sorted(snaplist)
+
+ return snaplist
+
+
+def compare_data(snaplist):
+ """
+ Create and save the plot
+ """
+
+ HI = []
+ HII = []
+ HeI = []
+ HeII = []
+ HeIII = []
+
+ for filename in snaplist:
+ data = swiftsimio.load(filename)
+
+ mXHI = data.gas.ion_mass_fractions.HI * data.gas.masses
+ mXHII = data.gas.ion_mass_fractions.HII * data.gas.masses
+ mXHeI = data.gas.ion_mass_fractions.HeI * data.gas.masses
+ mXHeII = data.gas.ion_mass_fractions.HeII * data.gas.masses
+ mXHeIII = data.gas.ion_mass_fractions.HeIII * data.gas.masses
+
+ HI.append(mXHI.sum())
+ HII.append(mXHII.sum())
+ HeI.append(mXHeI.sum())
+ HeII.append(mXHeII.sum())
+ HeIII.append(mXHeIII.sum())
+
+ plt.figure()
+ plt.plot(range(len(snaplist)), HI, label="HI total mass")
+ plt.plot(range(len(snaplist)), HII, label="HII total mass")
+ plt.plot(range(len(snaplist)), HeI, label="HeI total mass")
+ plt.plot(range(len(snaplist)), HeII, label="HeII total mass")
+ plt.plot(range(len(snaplist)), HeIII, label="HeIII total mass")
+ plt.legend()
+
+ # plt.show()
+ plt.tight_layout()
+ plt.savefig("total_abundancies.png", dpi=200)
+
+ return
+
+
+if __name__ == "__main__":
+
+ snaplist = get_snapshot_list(snapshot_base)
+ compare_data(snaplist)
diff --git a/examples/RadiativeTransferTests/IonizationEquilibriumICSetupTest/makeIC.py b/examples/RadiativeTransferTests/IonizationEquilibriumICSetupTest/makeIC.py
index a57098b4d22284d9001baa83723af170d395ec1d..4506f86e26970047a6b2a69c0c9701d9c66d2804 100755
--- a/examples/RadiativeTransferTests/IonizationEquilibriumICSetupTest/makeIC.py
+++ b/examples/RadiativeTransferTests/IonizationEquilibriumICSetupTest/makeIC.py
@@ -11,7 +11,6 @@ from swiftsimio import Writer
import unyt
import numpy as np
import h5py
-from matplotlib import pyplot as plt
# define unit system to use
unitsystem = unyt.unit_systems.cgs_unit_system
@@ -33,7 +32,7 @@ n_p = 10000
outputfilename = "ionization_equilibrium_test.hdf5"
# particle positions
-xp = unyt.unyt_array(np.zeros((n_p, 3), dtype=np.float), boxsize.units)
+xp = unyt.unyt_array(np.zeros((n_p, 3), dtype=np.float64), boxsize.units)
dx = boxsize / n_p
for i in range(n_p):
xp[i, 0] = (i + 0.5) * dx
@@ -42,7 +41,7 @@ w = Writer(unyt.unit_systems.cgs_unit_system, boxsize, dimension=1)
w.gas.coordinates = xp
w.gas.velocities = np.zeros(xp.shape) * (unyt.cm / unyt.s)
-w.gas.masses = np.ones(xp.shape[0], dtype=np.float) * 1000 * unyt.g
+w.gas.masses = np.ones(xp.shape[0], dtype=np.float64) * 1000 * unyt.g
w.gas.internal_energy = (
np.logspace(np.log10(umin.v), np.log10(umax.v), n_p) * unyt.erg / unyt.g
)
diff --git a/examples/RadiativeTransferTests/IonizationEquilibriumICSetupTest/plotSolution.py b/examples/RadiativeTransferTests/IonizationEquilibriumICSetupTest/plotSolution.py
index 39400b07008e2aa87bb8e78cc5dbb30897b1137f..d744fc0299f589e41cda3b87b36c589df46dde82 100755
--- a/examples/RadiativeTransferTests/IonizationEquilibriumICSetupTest/plotSolution.py
+++ b/examples/RadiativeTransferTests/IonizationEquilibriumICSetupTest/plotSolution.py
@@ -70,7 +70,7 @@ Tmax = T.v.max()
u_theory, mu_theory, T_theory, XHI_theory, XHII_theory, XHeI_theory, XHeII_theory, XHeIII_theory = np.loadtxt(
- "IonizationEquilibriumICSetupTestReference.txt", dtype=float, unpack=True
+ "IonizationEquilibriumICSetupTestReference.txt", dtype=np.float64, unpack=True
)
# ------------------
diff --git a/examples/RadiativeTransferTests/RandomizedBox_3D/README b/examples/RadiativeTransferTests/RandomizedBox_3D/README
index ef80ebb3314670dbc6692e92fa6c042fb298f69f..03d686ed719ea8693409699cae4f505362f04a73 100644
--- a/examples/RadiativeTransferTests/RandomizedBox_3D/README
+++ b/examples/RadiativeTransferTests/RandomizedBox_3D/README
@@ -1,6 +1,6 @@
This test case contains only gas and stars in a periodic box. The contents have
no physical meaning, and are background particles on a regular grid overlapping
-a randomly sampled sine wave particle distribution. The purpose of this test are
+a randomly sampled sine wave particle distribution. The purposes of this test are
as follows:
- Neither the stars nor the gas particles are updated all at the same time.
@@ -12,3 +12,6 @@ as follows:
You can and should run the same sanity check scripts as in ../UniformBox_3D to
check that everything's right.
+
+Compile swift with:
+ --with-rt=GEAR_1 --with-rt-riemann-solver=GLF --with-hydro-dimension=3 --with-hydro=gizmo-mfv --with-riemann-solver=hllc --with-stars=GEAR --with-feedback=none --enable-debugging-checks --with-grackle=$GRACKLE_ROOT
diff --git a/examples/RadiativeTransferTests/RandomizedBox_3D/makeIC.py b/examples/RadiativeTransferTests/RandomizedBox_3D/makeIC.py
index 452390e38ea28016421a7e43969c13bb180630b7..3b1850a6dcb092a276d73684e939076e1bbe251d 100755
--- a/examples/RadiativeTransferTests/RandomizedBox_3D/makeIC.py
+++ b/examples/RadiativeTransferTests/RandomizedBox_3D/makeIC.py
@@ -31,8 +31,6 @@ from swiftsimio.units import cosmo_units
import unyt
import numpy as np
-from scipy import special as sp
-from matplotlib import pyplot as plt
np.random.seed(666)
@@ -60,9 +58,9 @@ for i in range(n_p):
y = (j + 0.501) * dx
for k in range(n_p):
z = (k + 0.501) * dx
- xp.append(np.array([x, y, z], dtype=np.float))
+ xp.append(np.array([x, y, z], dtype=np.float64))
xp = np.array(xp)
-velp = np.zeros((n_p ** 3, 3), dtype=np.float)
+velp = np.zeros((n_p ** 3, 3), dtype=np.float64)
for i in range(n_s):
x = (i + 0.001) * ds
@@ -70,9 +68,9 @@ for i in range(n_s):
y = (j + 0.001) * ds
for k in range(n_s):
z = (k + 0.001) * ds
- xs.append(np.array([x, y, z], dtype=np.float))
+ xs.append(np.array([x, y, z], dtype=np.float64))
xs = np.array(xs)
-vels = np.zeros((n_s ** 3, 3), dtype=np.float)
+vels = np.zeros((n_s ** 3, 3), dtype=np.float64)
amplitude = 0.5
@@ -87,9 +85,9 @@ def sine(x, amplitude=1.0):
def sample(n):
samples = 0
- keep = np.zeros((n, 3), dtype=np.float)
+ keep = np.zeros((n, 3), dtype=np.float64)
while samples < n:
- sample = np.zeros(3, dtype=np.float)
+ sample = np.zeros(3, dtype=np.float64)
found = False
while not found:
@@ -155,11 +153,13 @@ w.gas.coordinates = xp
w.stars.coordinates = xs
w.gas.velocities = vp
w.stars.velocities = vs
-w.gas.masses = np.ones(xp.shape[0], dtype=np.float) * 1e6 * unyt.msun
+w.gas.masses = np.ones(xp.shape[0], dtype=np.float64) * 1e6 * unyt.msun
w.stars.masses = np.random.uniform(1e8, 1e10, size=xs.shape[0]) * unyt.msun
# Generate internal energy corresponding to 10^4 K
w.gas.internal_energy = (
- np.ones(xp.shape[0], dtype=float) * (1e4 * unyt.kb * unyt.K) / (1e6 * unyt.msun)
+ np.ones(xp.shape[0], dtype=np.float64)
+ * (1e4 * unyt.kb * unyt.K)
+ / (1e6 * unyt.msun)
)
diff --git a/examples/RadiativeTransferTests/RandomizedBox_3D/plotRadiationProjection.py b/examples/RadiativeTransferTests/RandomizedBox_3D/plotRadiationProjection.py
index e12c63ac7b231826dfe63e186c6076e768abb5de..f6edafe666937837d93eebab6b250a8798aa359e 100755
--- a/examples/RadiativeTransferTests/RandomizedBox_3D/plotRadiationProjection.py
+++ b/examples/RadiativeTransferTests/RandomizedBox_3D/plotRadiationProjection.py
@@ -17,9 +17,10 @@ import unyt
from matplotlib import pyplot as plt
import matplotlib as mpl
from mpl_toolkits.axes_grid1 import make_axes_locatable
-from matplotlib.colors import SymLogNorm, LogNorm
+from matplotlib.colors import SymLogNorm
# Parameters users should/may tweak
+# -----------------------------------
# plot all groups and all photon quantities
plot_all_data = True
@@ -158,9 +159,6 @@ def plot_photons(filename, energy_boundaries=None, flux_boundaries=None):
ax = fig.add_subplot(ngroups, 4, g * 4 + 1)
if energy_boundaries is not None:
- # imshow_kwargs["norm"] = None
- # imshow_kwargs["vmin"] = energy_boundaries[g][0]
- # imshow_kwargs["vmax"] = energy_boundaries[g][1]
imshow_kwargs["norm"] = SymLogNorm(
vmin=energy_boundaries[g][0],
vmax=energy_boundaries[g][1],
@@ -259,12 +257,6 @@ def plot_photons(filename, energy_boundaries=None, flux_boundaries=None):
if energy_boundaries is not None:
imshow_kwargs["vmin"] = energy_boundaries[g][0]
imshow_kwargs["vmax"] = energy_boundaries[g][1]
- # imshow_kwargs["norm"] = SymLogNorm(
- # vmin=energy_boundaries[g][0],
- # vmax=energy_boundaries[g][1],
- # linthresh = 1e-2,
- # base=10
- # )
im = ax.imshow(photon_map.T, **imshow_kwargs)
set_colorbar(ax, im)
ax.set_title("Group {0:2d}".format(g + 1))
@@ -331,7 +323,6 @@ def get_minmax_vals(snaplist):
dirmin = []
dirmax = []
for direction in ["X", "Y", "Z"]:
- new_attribute_str = "radiation_flux" + str(g + 1) + direction
f = getattr(data.gas.photon_fluxes, "Group" + str(g + 1) + direction)
dirmin.append(f.min())
dirmax.append(f.max())
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_2D/README b/examples/RadiativeTransferTests/StromgrenSphere_2D/README
index 0a0f3010ca265be9e4433fe7fcbc3762425c60bf..e9d6f81934406afab8a03d50f432e83736e5ed1a 100644
--- a/examples/RadiativeTransferTests/StromgrenSphere_2D/README
+++ b/examples/RadiativeTransferTests/StromgrenSphere_2D/README
@@ -1,22 +1,25 @@
Strömgen Sphere example in 2D
+-----------------------------
-For now, there are no gas interactions, and we use this
-example to perform tests and checks on photon injection
-and their propagation. The example contains 1 star in a
-glass hydro distribution in 2D. (State 09/2021)
+This directory contains two examples in one:
-Intented to compile with :
- --with-rt=GEAR_1 --with-hydro-dimension=2 --with-rt-riemann-solver=GLF
+ - run a Strömgren sphere example, where a single central
+ star ionizes the surrounding medium in a 2D plane.
+ To run this example, use the provided `run.sh` script.
+ This script will then make use of the `makeIC.py` and
+ `plotSolution.py` script.
-Scripts in this directory:
+ - run a propagation test of photons emitted from a single
+ central source in an otherwise uniform 2D plane.
+ To run this example, use the provided `runPropagationTest.sh` script.
+ This script will then make use of the `makePropagationTestIC.py` and
+ `plotPhotonPropagationCheck.py` script.
- - plotRadiationProjection.py:
- plot photon energies and fluxes as a 2D image
+Additional scripts:
+ - `plotRadiationProjection.py`: Plots a projection of the radiation
+ quantities (energies, fluxes). NOTE: you might need to change the
+ 'snapshot_base' variable at the top of the script depending on which
+ solutions you want to plot.
- - plotPhotonPropagationCheck.py:
- Plots radiation energies and fluxes as function of
- radius centered on the injecting star.
- NOTE: Make sure that you use a photon group that
- doesn't interact with the gas to check the
- propagation properly. You'll have to hard-code
- the index of the photon group in the script itself.
+To use the GEAR RT model, compile with :
+ --with-rt=GEAR_1 --with-rt-riemann-solver=GLF --with-hydro-dimension=2 --with-hydro=gizmo-mfv --with-riemann-solver=hllc --with-stars=GEAR --with-feedback=none --with-grackle=$GRACKLE_ROOT
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_2D/makeIC.py b/examples/RadiativeTransferTests/StromgrenSphere_2D/makeIC.py
index f57ee898b7fe7427734699d33cfbb67188e60e0c..3e39b9248a6e2fef65db51cbebd9d9d8fb0fa0ca 100755
--- a/examples/RadiativeTransferTests/StromgrenSphere_2D/makeIC.py
+++ b/examples/RadiativeTransferTests/StromgrenSphere_2D/makeIC.py
@@ -130,15 +130,8 @@ def get_number_densities_array(Temp, XH, XHe):
# n_He = X_He * rho_gas / m_He = (1 - X_H) / (4 X_H) * n_H
# = X_He / 4(1 - X_He) * nH = y * nH
- # if XH == 0:
- # nH = 0.0
- # nHe = 1.0
- # else:
- # nH = XH
- # nHe = XHe / 4
-
- nH = np.zeros(XH.shape, dtype=float)
- nHe = np.zeros(XH.shape, dtype=float)
+ nH = np.zeros(XH.shape, dtype=np.float64)
+ nHe = np.zeros(XH.shape, dtype=np.float64)
mask = XH == 0
nH[mask] = 0.0
@@ -148,13 +141,13 @@ def get_number_densities_array(Temp, XH, XHe):
nH[inv_mask] = XH[inv_mask]
nHe[inv_mask] = 0.25 * XHe[inv_mask]
- # NOTE: This is not the ionization threshold!
- nH0 = np.zeros(XH.shape, dtype=float)
- nHp = np.zeros(XH.shape, dtype=float)
- nHe0 = np.zeros(XH.shape, dtype=float)
- nHep = np.zeros(XH.shape, dtype=float)
- nHepp = np.zeros(XH.shape, dtype=float)
+ nH0 = np.zeros(XH.shape, dtype=np.float64)
+ nHp = np.zeros(XH.shape, dtype=np.float64)
+ nHe0 = np.zeros(XH.shape, dtype=np.float64)
+ nHep = np.zeros(XH.shape, dtype=np.float64)
+ nHepp = np.zeros(XH.shape, dtype=np.float64)
+ # NOTE: This is not the ionization threshold!
neutral = Temp < 5000 * unyt.K
nH0[neutral] = nH[neutral]
@@ -291,8 +284,8 @@ def mean_molecular_weight(XH0, XHp, XHe0, XHep, XHepp):
if __name__ == "__main__":
- glass = h5py.File("glassPlane_64.hdf5", "r")
- # glass = h5py.File("glassPlane_128.hdf5", "r")
+ # glass = h5py.File("glassPlane_64.hdf5", "r")
+ glass = h5py.File("glassPlane_128.hdf5", "r")
parts = glass["PartType0"]
xp = parts["Coordinates"][:]
h = parts["SmoothingLength"][:]
@@ -340,8 +333,8 @@ if __name__ == "__main__":
Mtot = rhoH * edgelen ** 3
mpart = Mtot / xp.shape[0]
mpart = mpart.to(cosmo_units["mass"])
- w.gas.masses = np.ones(xp.shape[0], dtype=np.float) * mpart
- w.stars.masses = np.ones(xs.shape[0], dtype=np.float) * mpart
+ w.gas.masses = np.ones(xp.shape[0], dtype=np.float64) * mpart
+ w.stars.masses = np.ones(xs.shape[0], dtype=np.float64) * mpart
# get gas internal energy for a given temperature and composition
XH = 1.0 # hydrogen mass fraction
@@ -351,6 +344,6 @@ if __name__ == "__main__":
mu = mean_molecular_weight(XHI, XHII, XHeI, XHeII, XHeIII)
internal_energy = internal_energy(T, mu)
- w.gas.internal_energy = np.ones(xp.shape[0], dtype=np.float) * internal_energy
+ w.gas.internal_energy = np.ones(xp.shape[0], dtype=np.float64) * internal_energy
w.write("stromgrenSphere-2D.hdf5")
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_2D/makePropagationTestIC.py b/examples/RadiativeTransferTests/StromgrenSphere_2D/makePropagationTestIC.py
new file mode 100755
index 0000000000000000000000000000000000000000..f1bc41207c606706d9e1abf80802b54b2e6b6c0c
--- /dev/null
+++ b/examples/RadiativeTransferTests/StromgrenSphere_2D/makePropagationTestIC.py
@@ -0,0 +1,57 @@
+#!/usr/bin/env python3
+
+# ---------------------------------------------------------------------
+# Add a single star in the center of a glass distribution
+# Intended for the photon propagation test.
+# ---------------------------------------------------------------------
+
+from swiftsimio import Writer
+import unyt
+import numpy as np
+import h5py
+
+
+glass = h5py.File("glassPlane_128.hdf5", "r")
+parts = glass["PartType0"]
+xp = parts["Coordinates"][:]
+h = parts["SmoothingLength"][:]
+glass.close()
+
+# replace the particle closest to the center
+# by the star
+r = np.sqrt(np.sum((0.5 - xp) ** 2, axis=1))
+rmin = np.argmin(r)
+xs = xp[rmin]
+xp = np.delete(xp, rmin, axis=0)
+h = np.delete(h, rmin)
+
+
+unitL = unyt.cm
+t_end = 1e-3 * unyt.s
+edgelen = unyt.c.to("cm/s") * t_end * 2.0
+edgelen = edgelen.to(unitL)
+boxsize = unyt.unyt_array([edgelen.v, edgelen.v, 0.0], unitL)
+
+xs = unyt.unyt_array(
+ [np.array([xs[0] * edgelen, xs[1] * edgelen, 0.0 * edgelen])], unitL
+)
+xp *= edgelen
+h *= edgelen
+
+
+w = Writer(unit_system=unyt.unit_systems.cgs_unit_system, box_size=boxsize, dimension=2)
+
+w.gas.coordinates = xp
+w.stars.coordinates = xs
+w.gas.velocities = np.zeros(xp.shape) * (unyt.cm / unyt.s)
+w.stars.velocities = np.zeros(xs.shape) * (unyt.cm / unyt.s)
+w.gas.masses = np.ones(xp.shape[0], dtype=np.float64) * 100 * unyt.g
+w.stars.masses = np.ones(xs.shape[0], dtype=np.float64) * 100 * unyt.g
+w.gas.internal_energy = (
+ np.ones(xp.shape[0], dtype=np.float64) * (300.0 * unyt.kb * unyt.K) / (unyt.g)
+)
+
+w.gas.smoothing_length = h
+w.stars.smoothing_length = w.gas.smoothing_length[:1]
+
+w.write("propagationTest-2D.hdf5")
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_2D/plotPhotonPropagationCheck.py b/examples/RadiativeTransferTests/StromgrenSphere_2D/plotPhotonPropagationCheck.py
index c8eb62f5dbe47841e957924b8e9d7333f9f32c52..548feee03bd4e94ccc3ce089cf1e5125f2e5d988 100755
--- a/examples/RadiativeTransferTests/StromgrenSphere_2D/plotPhotonPropagationCheck.py
+++ b/examples/RadiativeTransferTests/StromgrenSphere_2D/plotPhotonPropagationCheck.py
@@ -16,21 +16,26 @@
# correctly.
# ----------------------------------------------------------------------
-import sys, os, gc
+import sys
+import os
+import gc
import swiftsimio
import numpy as np
import unyt
from matplotlib import pyplot as plt
import matplotlib as mpl
from scipy import stats
-from scipy.optimize import curve_fit
# Parameters users should/may tweak
# snapshot basename
-snapshot_base = "output"
+snapshot_base = "propagation_test"
+
+# additional anisotropy estimate plot?
+plot_anisotropy_estimate = False
# which photon group to use.
+# NOTE: array index, not group number (which starts at 1 for GEAR)
group_index = 0
scatterplot_kwargs = {
@@ -55,7 +60,7 @@ except IndexError:
mpl.rcParams["text.usetex"] = True
-def analytical_intgrated_energy_solution(L, time, r, rmax):
+def analytical_integrated_energy_solution(L, time, r, rmax):
"""
Compute analytical solution for the sum of the energy
in bins for given injection rate <L> at time <time>
@@ -169,6 +174,7 @@ def plot_photons(filename, emin, emax, fmin, fmax):
r_expect = meta.time * meta.reduced_lightspeed
use_const_emission_rates = bool(meta.parameters["GEARRT:use_const_emission_rates"])
+ L = None
if use_const_emission_rates:
# read emission rate parameter as string
emissionstr = meta.parameters["GEARRT:star_emission_rates_LSol"].decode("utf-8")
@@ -186,7 +192,11 @@ def plot_photons(filename, emin, emax, fmin, fmax):
const_emission_rates = unyt.unyt_array(emlist, unyt.L_Sun)
L = const_emission_rates[group_index]
- fig = plt.figure(figsize=(5 * 4, 5.5), dpi=200)
+ if plot_anisotropy_estimate:
+ ncols = 4
+ else:
+ ncols = 3
+ fig = plt.figure(figsize=(5 * ncols, 5.5), dpi=200)
nbins = 100
r_bin_edges = np.linspace(0.5 * edgelen * 1e-2, 0.507 * edgelen, nbins + 1)
@@ -204,8 +214,6 @@ def plot_photons(filename, emin, emax, fmin, fmax):
Fy = getattr(data.gas.photon_fluxes, "Group" + str(group_index + 1) + "Y")
fmag = np.sqrt(Fx ** 2 + Fy ** 2)
- sum_fmag = fmag.sum()
- max_fmag = fmag.max()
particle_count, _ = np.histogram(
r,
bins=r_analytical_bin_edges,
@@ -220,7 +228,7 @@ def plot_photons(filename, emin, emax, fmin, fmax):
# ------------------------
# Plot photon energies
# ------------------------
- ax1 = fig.add_subplot(1, 4, 1)
+ ax1 = fig.add_subplot(1, ncols, 1)
ax1.set_title("Particle Radiation Energies")
ax1.set_ylabel("Photon Energy [$" + energy_units_str + "$]")
@@ -278,7 +286,7 @@ def plot_photons(filename, emin, emax, fmin, fmax):
# ------------------------------
# Plot binned photon energies
# ------------------------------
- ax2 = fig.add_subplot(1, 4, 2)
+ ax2 = fig.add_subplot(1, ncols, 2)
ax2.set_title("Total Radiation Energy in radial bins")
ax2.set_ylabel("Total Photon Energy [$" + energy_units_str + "$]")
@@ -301,7 +309,7 @@ def plot_photons(filename, emin, emax, fmin, fmax):
if use_const_emission_rates:
# plot entire expected solution
# Note: you need to use the same bins as for the actual results
- rA, EA = analytical_intgrated_energy_solution(L, time, r_bin_edges, r_expect)
+ rA, EA = analytical_integrated_energy_solution(L, time, r_bin_edges, r_expect)
ax2.plot(
rA,
@@ -323,7 +331,7 @@ def plot_photons(filename, emin, emax, fmin, fmax):
# ------------------------------
# Plot photon fluxes
# ------------------------------
- ax3 = fig.add_subplot(1, 4, 3)
+ ax3 = fig.add_subplot(1, ncols, 3)
ax3.set_title("Particle Radiation Flux Magnitudes")
ax3.set_ylabel("Photon Flux Magnitude [$" + flux_units_str + "$]")
@@ -387,55 +395,57 @@ def plot_photons(filename, emin, emax, fmin, fmax):
# Plot photon flux sum
# ------------------------------
- ax4 = fig.add_subplot(1, 4, 4)
- ax4.set_title("Vectorial Sum of Radiation Flux in radial bins")
- ax4.set_ylabel("[1]")
+ if plot_anisotropy_estimate:
- fmag_sum_bin, _, _ = stats.binned_statistic(
- r,
- fmag,
- statistic="sum",
- bins=r_bin_edges,
- range=(r_bin_edges[0], r_bin_edges[-1]),
- )
- mask_sum = fmag_sum_bin > 0
- fmag_max_bin, _, _ = stats.binned_statistic(
- r,
- fmag,
- statistic="max",
- bins=r_bin_edges,
- range=(r_bin_edges[0], r_bin_edges[-1]),
- )
- mask_max = fmag_max_bin > 0
- Fx_sum_bin, _, _ = stats.binned_statistic(
- r,
- Fx,
- statistic="sum",
- bins=r_bin_edges,
- range=(r_bin_edges[0], r_bin_edges[-1]),
- )
- Fy_sum_bin, _, _ = stats.binned_statistic(
- r,
- Fy,
- statistic="sum",
- bins=r_bin_edges,
- range=(r_bin_edges[0], r_bin_edges[-1]),
- )
- F_sum_bin = np.sqrt(Fx_sum_bin ** 2 + Fy_sum_bin ** 2)
+ ax4 = fig.add_subplot(1, ncols, 4)
+ ax4.set_title("Vectorial Sum of Radiation Flux in radial bins")
+ ax4.set_ylabel("[1]")
- ax4.plot(
- r_bin_centres[mask_sum],
- F_sum_bin[mask_sum] / fmag_sum_bin[mask_sum],
- **lineplot_kwargs,
- label="$\left| \sum_{i \in \mathrm{particles \ in \ bin}} \mathbf{F}_i \\right| $ / $\sum_{i \in \mathrm{particles \ in \ bin}} \left| \mathbf{F}_{i} \\right| $",
- )
- ax4.plot(
- r_bin_centres[mask_max],
- F_sum_bin[mask_max] / fmag_max_bin[mask_max],
- **lineplot_kwargs,
- linestyle="--",
- label="$\left| \sum_{i \in \mathrm{particles \ in \ bin}} \mathbf{F}_i \\right| $ / $\max_{i \in \mathrm{particles \ in \ bin}} \left| \mathbf{F}_{i} \\right| $",
- )
+ fmag_sum_bin, _, _ = stats.binned_statistic(
+ r,
+ fmag,
+ statistic="sum",
+ bins=r_bin_edges,
+ range=(r_bin_edges[0], r_bin_edges[-1]),
+ )
+ mask_sum = fmag_sum_bin > 0
+ fmag_max_bin, _, _ = stats.binned_statistic(
+ r,
+ fmag,
+ statistic="max",
+ bins=r_bin_edges,
+ range=(r_bin_edges[0], r_bin_edges[-1]),
+ )
+ mask_max = fmag_max_bin > 0
+ Fx_sum_bin, _, _ = stats.binned_statistic(
+ r,
+ Fx,
+ statistic="sum",
+ bins=r_bin_edges,
+ range=(r_bin_edges[0], r_bin_edges[-1]),
+ )
+ Fy_sum_bin, _, _ = stats.binned_statistic(
+ r,
+ Fy,
+ statistic="sum",
+ bins=r_bin_edges,
+ range=(r_bin_edges[0], r_bin_edges[-1]),
+ )
+ F_sum_bin = np.sqrt(Fx_sum_bin ** 2 + Fy_sum_bin ** 2)
+
+ ax4.plot(
+ r_bin_centres[mask_sum],
+ F_sum_bin[mask_sum] / fmag_sum_bin[mask_sum],
+ **lineplot_kwargs,
+ label="$\left| \sum_{i \in \mathrm{particles \ in \ bin}} \mathbf{F}_i \\right| $ / $\sum_{i \in \mathrm{particles \ in \ bin}} \left| \mathbf{F}_{i} \\right| $",
+ )
+ ax4.plot(
+ r_bin_centres[mask_max],
+ F_sum_bin[mask_max] / fmag_max_bin[mask_max],
+ **lineplot_kwargs,
+ linestyle="--",
+ label="$\left| \sum_{i \in \mathrm{particles \ in \ bin}} \mathbf{F}_i \\right| $ / $\max_{i \in \mathrm{particles \ in \ bin}} \left| \mathbf{F}_{i} \\right| $",
+ )
# -------------------------------------------
# Cosmetics that all axes have in common
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_2D/plotRadiationProjection.py b/examples/RadiativeTransferTests/StromgrenSphere_2D/plotRadiationProjection.py
index 148c5e0c032e3ead59689881a3ec9ed1eca2310d..ef6298538bf3f9abd840a49a2bf341d537071776 100755
--- a/examples/RadiativeTransferTests/StromgrenSphere_2D/plotRadiationProjection.py
+++ b/examples/RadiativeTransferTests/StromgrenSphere_2D/plotRadiationProjection.py
@@ -11,7 +11,6 @@
import sys
import os
import swiftsimio
-import numpy as np
import gc
import unyt
from matplotlib import pyplot as plt
@@ -23,8 +22,10 @@ from matplotlib.colors import SymLogNorm
# plot all groups and all photon quantities
plot_all_data = True
-# snapshot basename
-snapshot_base = "output"
+
+# plotting propagation tests or Stromgren sphere?
+do_stromgren_sphere = True
+
# fancy up the plots a bit?
fancy = True
@@ -34,12 +35,24 @@ imshow_kwargs = {"origin": "lower", "cmap": "viridis"}
# parameters for swiftsimio projections
projection_kwargs = {"resolution": 1024, "parallel": True}
+# snapshot basename
+snapshot_base = "propagation_test"
+
# Set Units of your choice
energy_units = unyt.erg
energy_units_str = "\\rm{erg}"
flux_units = 1e10 * energy_units / unyt.cm ** 2 / unyt.s
flux_units_str = "10^{10} \\rm{erg} \\ \\rm{cm}^{-2} \\ \\rm{s}^{-1}"
time_units = unyt.s
+
+if do_stromgren_sphere:
+ snapshot_base = "output"
+
+ energy_units = 1e50 * unyt.erg
+ energy_units_str = "10^{50} \\rm{erg}"
+ flux_units = 1e50 * unyt.erg / unyt.kpc ** 2 / unyt.Gyr
+ flux_units_str = "10^{60} \\rm{erg} \\ \\rm{kpc}^{-2} \\ \\rm{Gyr}^{-1}"
+ time_units = unyt.Myr
# -----------------------------------------------------------------------
@@ -76,6 +89,10 @@ def get_snapshot_list(snapshot_basename="output"):
quit(1)
snaplist.append(fname)
+ if len(snaplist) == 0:
+ print("Didn't find any snapshots with basename '", snapshot_basename, "'")
+ quit(1)
+
return snaplist
@@ -139,7 +156,7 @@ def plot_photons(filename, energy_boundaries=None, flux_boundaries=None):
if plot_all_data:
fig = plt.figure(figsize=(5 * 3, 5.05 * ngroups), dpi=200)
- figname = filename[:-5] + "-all-quantities.png"
+ figname = filename[:-5] + "-radiation.png"
for g in range(ngroups):
@@ -296,7 +313,6 @@ def get_minmax_vals(snaplist):
dirmin = []
dirmax = []
for direction in ["X", "Y"]:
- new_attribute_str = "radiation_flux" + str(g + 1) + direction
f = getattr(data.gas.photon_fluxes, "Group" + str(g + 1) + direction)
dirmin.append(f.min())
dirmax.append(f.max())
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_2D/plotSolution.py b/examples/RadiativeTransferTests/StromgrenSphere_2D/plotSolution.py
index 9417b943d80d7f5871d38893647b415494a6a1be..06257017567beb26027544d9c9ce09689e482f7f 100755
--- a/examples/RadiativeTransferTests/StromgrenSphere_2D/plotSolution.py
+++ b/examples/RadiativeTransferTests/StromgrenSphere_2D/plotSolution.py
@@ -8,14 +8,12 @@
import sys
import os
import swiftsimio
-import numpy as np
import gc
import unyt
from matplotlib import pyplot as plt
import matplotlib as mpl
from mpl_toolkits.axes_grid1 import make_axes_locatable
-from matplotlib.colors import SymLogNorm, LogNorm
-from swiftsimio.visualisation.slice import slice_gas
+from matplotlib.colors import LogNorm
# Parameters users should/may tweak
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_2D/propagationTest-2D.yml b/examples/RadiativeTransferTests/StromgrenSphere_2D/propagationTest-2D.yml
new file mode 100644
index 0000000000000000000000000000000000000000..3b0e7190c40263969969702c4530900002274da4
--- /dev/null
+++ b/examples/RadiativeTransferTests/StromgrenSphere_2D/propagationTest-2D.yml
@@ -0,0 +1,60 @@
+MetaData:
+ run_name: StromgrenSpherePropagationTest2D
+
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.
+ UnitLength_in_cgs: 1.
+ UnitVelocity_in_cgs: 1.
+ UnitCurrent_in_cgs: 1.
+ UnitTemp_in_cgs: 1.
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 0.001 # 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: 2.e-04 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: propagation_test # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 0.0001 # Time between snapshots
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ time_first: 0.
+ delta_time: 5e-4 # 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.6 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 10. # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./propagationTest-2D.hdf5 # The file to read
+ periodic: 1 # peridioc ICs. Keep them periodic so we don't loose photon energy. TODO: CHANGE LATER WHEN YOU ACTUALLY DO GAS INTERACTIONS
+
+Scheduler:
+ max_top_level_cells: 12
+
+GEARRT:
+ f_reduce_c: 1. # reduce the speed of light for the RT solver by multiplying c with this factor
+ CFL_condition: 0.99
+ photon_groups_Hz: [] # Photon frequency group bin edges in Hz. Needs to be 1 less than the number of groups (N) requested during the configuration (--with-RT=GEAR_N).
+ use_const_emission_rates: 1 # (Optional) use constant emission rates for stars as defined with star_emission_rates_erg_LSol parameter
+ star_emission_rates_LSol: [1e-28] # (Optional) constant star emission rates for each photon frequency group to use if use_constant_emission_rates is set.
+ hydrogen_mass_fraction: 0.76 # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
+ stellar_spectrum_type: 0 # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
+ stellar_spectrum_const_max_frequency_Hz: 1.e17 # (Conditional) if stellar_spectrum_type=0, use this maximal frequency for the constant spectrum.
+ set_initial_ionization_mass_fractions: 1 # (Optional) manually overwrite initial mass fraction of each species (using the values you set below)
+ mass_fraction_HI: 0.76 # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HI mass fractions to this value
+ mass_fraction_HII: 0. # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HII mass fractions to this value
+ mass_fraction_HeI: 0.24 # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeI mass fractions to this value
+ mass_fraction_HeII: 0. # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeII mass fractions to this value
+ mass_fraction_HeIII: 0. # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeIII mass fractions to this value
+ skip_thermochemistry: 1 # skip thermochemistry.
+ stars_max_timestep: 7.812500e-06 # update stars every step!
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_2D/run.sh b/examples/RadiativeTransferTests/StromgrenSphere_2D/run.sh
index 4db90a0b079f3c92c4f5bdb079b13edb709707c5..74153c46bdb51790cbca1dfa78f3832031ca5c63 100755
--- a/examples/RadiativeTransferTests/StromgrenSphere_2D/run.sh
+++ b/examples/RadiativeTransferTests/StromgrenSphere_2D/run.sh
@@ -1,5 +1,9 @@
#!/bin/bash
+#---------------------------------------
+# Runs the Stromgren Sphere example
+#---------------------------------------
+
# make run.sh fail if a subcommand fails
set -e
set -o pipefail
@@ -18,11 +22,12 @@ fi
# Run SWIFT with RT
../../swift \
--hydro --threads=4 --stars --external-gravity \
- --feedback --radiation --fpe \
+ --feedback --radiation \
+ --fpe \
stromgrenSphere-2D.yml 2>&1 | tee output.log
# Plot the photon propagation checks.
# Make sure you set the correct photon group to plot
# inside the script
-python3 ./plotPhotonPropagationCheck.py
+python3 ./plotSolution.py
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_2D/runPropagationTest.sh b/examples/RadiativeTransferTests/StromgrenSphere_2D/runPropagationTest.sh
new file mode 100755
index 0000000000000000000000000000000000000000..83ec3c33c1a98e8640820e4139109c6c725e4673
--- /dev/null
+++ b/examples/RadiativeTransferTests/StromgrenSphere_2D/runPropagationTest.sh
@@ -0,0 +1,32 @@
+#!/bin/bash
+
+#---------------------------------------
+# Runs the Propagation Test example
+#---------------------------------------
+
+
+# make run.sh fail if a subcommand fails
+set -e
+set -o pipefail
+
+if [ ! -e glassPlane_128.hdf5 ]
+then
+ echo "Fetching initial glass file for Strömgen Sphere 2D example ..."
+ ./getGlass.sh
+fi
+
+if [ ! -f 'propagationTest-2D.hdf5' ]; then
+ echo "Generating ICs"
+ python3 makePropagationTestIC.py
+fi
+
+# Run SWIFT with RT
+../../swift \
+ --hydro --threads=4 --stars --external-gravity \
+ --feedback --radiation --fpe \
+ ./propagationTest-2D.yml 2>&1 | tee output.log
+
+# Plot the photon propagation checks.
+# Make sure you set the correct photon group to plot
+# inside the script
+python3 ./plotPhotonPropagationCheck.py
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_2D/stromgrenSphere-2D.yml b/examples/RadiativeTransferTests/StromgrenSphere_2D/stromgrenSphere-2D.yml
index 8b6cc7688869c60416672708509ff1c91e753f18..5e8d44eb2b873476d331e2af6db159b4568ad6e9 100644
--- a/examples/RadiativeTransferTests/StromgrenSphere_2D/stromgrenSphere-2D.yml
+++ b/examples/RadiativeTransferTests/StromgrenSphere_2D/stromgrenSphere-2D.yml
@@ -8,33 +8,24 @@ InternalUnitSystem:
UnitVelocity_in_cgs: 1.e5 # km/s
UnitCurrent_in_cgs: 1.
UnitTemp_in_cgs: 1. # K
- # UnitMass_in_cgs: 7.98841586e+23 # 1 M_Sol
- # UnitLength_in_cgs: 7.08567758e12 # kpc in cm
- # UnitVelocity_in_cgs: 113.e7 # km/s
- # UnitCurrent_in_cgs: 1.
- # UnitTemp_in_cgs: 123. # K
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 0.512
- # time_end: 4.882813e-03
- # time_end: 5.960465e-06
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: 5.960465e-07
-
# Parameters governing the snapshots
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.004 # Time between snapshots
+ delta_time: 0.001 # Time between snapshots
# Parameters governing the conserved quantities statistics
Statistics:
time_first: 0.
- delta_time: 1e-4 # Time between statistics output
+ delta_time: 1e-3 # Time between statistics output
# Parameters for the hydrodynamics scheme
SPH:
@@ -45,7 +36,7 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
file_name: ./stromgrenSphere-2D.hdf5 # The file to read
- periodic: 0 # periodic ICs. Keep them periodic so we don't loose photon energy. TODO: CHANGE LATER WHEN YOU ACTUALLY DO GAS INTERACTIONS
+ periodic: 0 # periodic ICs. Keep them periodic so we don't loose photon energy.
Scheduler:
max_top_level_cells: 12
@@ -58,8 +49,6 @@ GEARRT:
CFL_condition: 0.9
photon_groups_Hz: [3.288e15, 5.945e15, 13.157e15] # Photon frequency group bin edges in Hz. Needs to be 1 less than the number of groups (N) requested during the configuration (--with-RT=GEAR_N). Outer edges of zero and infinity are assumed.
use_const_emission_rates: 1 # (Optional) use constant emission rates for stars as defined with star_emission_rates_erg_LSol parameter
- # star_emission_rates_LSol: [1000000., 1000000., 1000000., 1000000.] # (Optional) constant star emission rates for each photon frequency group to use if use_constant_emission_rates is set, in units of Solar Luminosity.
- # star_emission_rates_LSol: [61996.44159, 61996.44159] # (Optional) constant star emission rates for each photon frequency group to use if use_constant_emission_rates is set, in units of Solar Luminosity.
star_emission_rates_LSol: [247985.76636, 247985.76636] # (Optional) constant star emission rates for each photon frequency group to use if use_constant_emission_rates is set, in units of Solar Luminosity.
hydrogen_mass_fraction: 1.00 # 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.
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_3D/README b/examples/RadiativeTransferTests/StromgrenSphere_3D/README
index da3fd209a4945dd051a4e1a9f9d74e7ac84a0ab1..a123177f03b2c06af367c719be442646deb62fef 100644
--- a/examples/RadiativeTransferTests/StromgrenSphere_3D/README
+++ b/examples/RadiativeTransferTests/StromgrenSphere_3D/README
@@ -1,23 +1,27 @@
Strömgen Sphere example in 3D
+-----------------------------
-For now, there are no gas interactions, and we use this
-example to perform tests and checks on photon injection
-and their propagation. The example contains 1 star in a
-glass hydro distribution in 3D. (State 09/2021)
+This directory contains two examples in one:
-Intented to compile with :
- --with-rt=GEAR_1 --with-hydro-dimension=3
+ - run a Strömgren sphere example, where a single central
+ star ionizes the surrounding medium in a box.
+ To run this example, use the provided `run.sh` script.
+ This script will then make use of the `makeIC.py` and
+ `plotSolution.py` script.
-Scripts in this directory:
+ - run a propagation test of photons emitted from a single
+ central source in an otherwise uniform box.
+ To run this example, use the provided `runPropagationTest.sh` script.
+ This script will then make use of the `makePropagationTestIC.py` and
+ `plotPhotonPropagationCheck.py` script.
- - plotRadiationProjection.py:
- plot photon energies and fluxes as a 2D projected
- image
- - plotPhotonPropagationCheck.py:
- Plots radiation energies and fluxes as function of
- radius centered on the injecting star.
- NOTE: Make sure that you use a photon group that
- doesn't interact with the gas to check the
- propagation properly. You'll have to hard-code
- the index of the photon group in the script itself.
+Additional scripts:
+ - `plotRadiationProjection.py`: Plots a projection of the radiation
+ quantities (energies, fluxes). NOTE: you might need to change the
+ 'snapshot_base' variable at the top of the script depending on which
+ solutions you want to plot.
+
+To use the GEAR RT model, compile with :
+ --with-rt=GEAR_1 --with-rt-riemann-solver=GLF --with-hydro=gizmo-mfv --with-riemann-solver=hllc --with-stars=GEAR --with-feedback=none --with-grackle=$GRACKLE_ROOT
+
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_3D/makeIC.py b/examples/RadiativeTransferTests/StromgrenSphere_3D/makeIC.py
index 384178ba209617144391170a1a95bd7d0f619402..dc34328b232ec98ee89b2107add90832167ae36a 100755
--- a/examples/RadiativeTransferTests/StromgrenSphere_3D/makeIC.py
+++ b/examples/RadiativeTransferTests/StromgrenSphere_3D/makeIC.py
@@ -5,52 +5,344 @@
# ---------------------------------------------------------------------
from swiftsimio import Writer
+from swiftsimio.units import cosmo_units
import unyt
import numpy as np
import h5py
+gamma = 5.0 / 3.0
-glass = h5py.File("glassCube_64.hdf5", "r")
-parts = glass["PartType0"]
-xp = parts["Coordinates"][:]
-h = parts["SmoothingLength"][:]
-glass.close()
-# replace the particle closest to the center
-# by the star
-r = np.sqrt(np.sum((0.5 - xp) ** 2, axis=1))
-rmin = np.argmin(r)
-xs = xp[rmin]
-xp = np.delete(xp, rmin, axis=0)
-h = np.delete(h, rmin)
+def get_number_densities(Temp, XH, XHe):
+ """
+ Compute the number densities of all species at a given
+ temperature following Katz, Hernquist, and Weinberg 1996
+ Temp: temperature [unyt quantity]
+ XH: total mass fraction of all hydrogen species (HI + HII)
+ XHe: total mass fraction of all helium species (HeI + HeII + HeIII)
+ """
-unitL = unyt.cm
-t_end = 1e-3 * unyt.s
-edgelen = unyt.c.to("cm/s") * t_end * 2.0
-edgelen = edgelen.to(unitL)
-boxsize = unyt.unyt_array([edgelen.v, edgelen.v, edgelen.v], unitL)
+ # n_H = X_H * rho_gas / m_H =
+ # n_He = X_He * rho_gas / m_He = (1 - X_H) / (4 X_H) * n_H
+ # = X_He / 4(1 - X_He) * nH = y * nH
-xs = unyt.unyt_array(
- [np.array([xs[0] * edgelen, xs[1] * edgelen, xs[2] * edgelen])], unitL
-)
-xp *= edgelen
-h *= edgelen
+ if XH == 0:
+ nH = 0.0
+ nHe = 1.0
+ else:
+ nH = XH
+ nHe = XHe / 4
+ # NOTE: This is not the ionization threshold!
+ if Temp < 5000 * unyt.K:
+ nH0 = nH
+ nHp = 0.0
+ nHe0 = nHe
+ nHep = 0.0
+ nHepp = 0.0
-w = Writer(unit_system=unyt.unit_systems.cgs_unit_system, box_size=boxsize, dimension=3)
+ else:
-w.gas.coordinates = xp
-w.stars.coordinates = xs
-w.gas.velocities = np.zeros(xp.shape) * (unyt.cm / unyt.s)
-w.stars.velocities = np.zeros(xs.shape) * (unyt.cm / unyt.s)
-w.gas.masses = np.ones(xp.shape[0], dtype=np.float) * 100 * unyt.g
-w.stars.masses = np.ones(xs.shape[0], dtype=np.float) * 100 * unyt.g
-w.gas.internal_energy = (
- np.ones(xp.shape[0], dtype=np.float) * (300.0 * unyt.kb * unyt.K) / (unyt.g)
-)
+ Temp.convert_to_cgs()
+ T = Temp.v
+ # Recombination rate for H+ in units of cm^3 s^-1
+ A_Hp = (
+ 8.40e-11
+ / np.sqrt(T)
+ * (T * 1e-3) ** (-0.2)
+ * 1.0
+ / (1.0 + (T * 1e-6) ** 0.7)
+ )
-w.gas.smoothing_length = h
-w.stars.smoothing_length = w.gas.smoothing_length[:1]
+ # Dielectronic recombination rate for He+ in units of cm^3 s^-1
+ A_d = (
+ 1.9e-3
+ / T ** 1.5
+ * np.exp(-470000.0 / T)
+ * (1.0 + 0.3 * np.exp(-94000.0 / T))
+ )
+ # Recombination rate for He+ in units of cm^3 s^-1
+ A_Hep = 1.5e-10 / T ** 0.6353
+ # Recombination rate for He++ in units of cm^3 s^-1
+ A_Hepp = (
+ 3.36e-10
+ / np.sqrt(T)
+ * (T * 1e-3) ** (-0.2)
+ * 1.0
+ / (1.0 + (T * 1e-6) ** 0.7)
+ )
+ # collisional ionization rate for H0 in units of cm^3 s^-1
+ # G_H0 = 1.17e-10 * np.sqrt(T) * np.exp(-157809.1 / T) * 1. / (1. + np.sqrt(T*1e-5))
+ G_H0 = (
+ 5.85e-11
+ * np.sqrt(T)
+ * np.exp(-157809.1 / T)
+ * 1.0
+ / (1.0 + np.sqrt(T * 1e-5))
+ )
+ # collisional ionization rate for He0 in units of cm^3 s^-1
+ G_He0 = (
+ 2.38e-11
+ * np.sqrt(T)
+ * np.exp(-285335.4 / T)
+ * 1.0
+ / (1.0 + np.sqrt(T * 1e-5))
+ )
+ # collisional ionization rate for He+ in units of cm^3 s^-1
+ G_Hep = (
+ 5.68e-12
+ * np.sqrt(T)
+ * np.exp(-631515.0 / T)
+ * 1.0
+ / (1.0 + np.sqrt(T * 1e-5))
+ )
-w.write("stromgrenSphere-3D.hdf5")
+ # Katz et al. 1996 eq. 33 - 38
+ # Note: We assume all photoionization rates to be zero.
+ # Also, here we don't care about the actual number density, i.e.
+ # about the volume, since it'll cancel out later when we compute
+ # the mass fractions.
+
+ nH0 = nH * A_Hp / (A_Hp + G_H0)
+ nHp = nH - nH0
+ nHep = nHe / (1.0 + (A_Hep + A_d) / G_He0 + G_Hep / A_Hepp)
+ nHe0 = nHep * (A_Hep + A_d) / G_He0
+ nHepp = nHep * G_Hep / A_Hepp
+
+ # electron density
+ ne = nHp + nHep + 2.0 * nHepp
+
+ return nH0, nHp, nHe0, nHep, nHepp, ne
+
+
+def get_number_densities_array(Temp, XH, XHe):
+ """
+ Compute the number densities of all species at a given
+ temperature following Katz, Hernquist, and Weinberg 1996
+
+ Temp: temperature [unyt array]
+ XH: total mass fraction of all hydrogen species (HI + HII)
+ XHe: total mass fraction of all helium species (HeI + HeII + HeIII)
+ """
+
+ # n_H = X_H * rho_gas / m_H =
+ # n_He = X_He * rho_gas / m_He = (1 - X_H) / (4 X_H) * n_H
+ # = X_He / 4(1 - X_He) * nH = y * nH
+
+ nH = np.zeros(XH.shape, dtype=np.float64)
+ nHe = np.zeros(XH.shape, dtype=np.float64)
+
+ mask = XH == 0
+ nH[mask] = 0.0
+ nHe[mask] = 1.0
+
+ inv_mask = np.logical_not(mask)
+ nH[inv_mask] = XH[inv_mask]
+ nHe[inv_mask] = 0.25 * XHe[inv_mask]
+
+ nH0 = np.zeros(XH.shape, dtype=np.float64)
+ nHp = np.zeros(XH.shape, dtype=np.float64)
+ nHe0 = np.zeros(XH.shape, dtype=np.float64)
+ nHep = np.zeros(XH.shape, dtype=np.float64)
+ nHepp = np.zeros(XH.shape, dtype=np.float64)
+
+ # NOTE: This is not the ionization threshold!
+ neutral = Temp < 5000 * unyt.K
+
+ nH0[neutral] = nH[neutral]
+ nHp[neutral] = 0.0
+ nHe0[neutral] = nHe[neutral]
+ nHep[neutral] = 0.0
+ nHepp[neutral] = 0.0
+
+ Temp.convert_to_cgs()
+ T = Temp.v
+ ionized = np.logical_not(neutral)
+
+ # Recombination rate for H+ in units of cm^3 s^-1
+ A_Hp = (
+ 8.40e-11 / np.sqrt(T) * (T * 1e-3) ** (-0.2) * 1.0 / (1.0 + (T * 1e-6) ** 0.7)
+ )
+ # Dielectronic recombination rate for He+ in units of cm^3 s^-1
+ A_d = 1.9e-3 / T ** 1.5 * np.exp(-470000.0 / T) * (1.0 + 0.3 * np.exp(-94000.0 / T))
+ # Recombination rate for He+ in units of cm^3 s^-1
+ A_Hep = 1.5e-10 / T ** 0.6353
+ # Recombination rate for He++ in units of cm^3 s^-1
+ A_Hepp = (
+ 3.36e-10 / np.sqrt(T) * (T * 1e-3) ** (-0.2) * 1.0 / (1.0 + (T * 1e-6) ** 0.7)
+ )
+ # collisional ionization rate for H0 in units of cm^3 s^-1
+ G_H0 = (
+ 5.85e-11 * np.sqrt(T) * np.exp(-157809.1 / T) * 1.0 / (1.0 + np.sqrt(T * 1e-5))
+ )
+ # collisional ionization rate for He0 in units of cm^3 s^-1
+ G_He0 = (
+ 2.38e-11 * np.sqrt(T) * np.exp(-285335.4 / T) * 1.0 / (1.0 + np.sqrt(T * 1e-5))
+ )
+ # collisional ionization rate for He+ in units of cm^3 s^-1
+ G_Hep = (
+ 5.68e-12 * np.sqrt(T) * np.exp(-631515.0 / T) * 1.0 / (1.0 + np.sqrt(T * 1e-5))
+ )
+
+ # Katz et al. 1996 eq. 33 - 38
+ # Note: We assume all photoionization rates to be zero.
+ # Also, here we don't care about the actual number density, i.e.
+ # about the volume, since it'll cancel out later when we compute
+ # the mass fractions.
+
+ nH0[ionized] = nH[ionized] * A_Hp[ionized] / (A_Hp[ionized] + G_H0[ionized])
+ nHp[ionized] = nH[ionized] - nH0[ionized]
+ nHep[ionized] = nHe[ionized] / (
+ 1.0
+ + (A_Hep[ionized] + A_d[ionized]) / G_He0[ionized]
+ + G_Hep[ionized] / A_Hepp[ionized]
+ )
+ nHe0[ionized] = nHep[ionized] * (A_Hep[ionized] + A_d[ionized]) / G_He0[ionized]
+ nHepp[ionized] = nHep[ionized] * G_Hep[ionized] / A_Hepp[ionized]
+
+ # electron density
+ ne = nHp + nHep + 2.0 * nHepp
+
+ return nH0, nHp, nHe0, nHep, nHepp, ne
+
+
+def get_mass_fractions(T, XH, XHe):
+ """
+ Compute the mass fractions of all species at a
+ given temperature
+
+ T: temperature
+ XH: total mass fraction of all hydrogen species (HI + HII)
+ XHe: total mass fraction of all helium species (HeI + HeII + HeIII)
+ """
+
+ # first get number densities
+ if isinstance(XH, np.ndarray):
+ nH0, nHp, nHe0, nHep, nHepp, ne = get_number_densities_array(T, XH, XHe)
+ else:
+ nH0, nHp, nHe0, nHep, nHepp, ne = get_number_densities(T, XH, XHe)
+
+ # now get mass denities in units of atomic mass units
+ mH0 = nH0
+ mHp = nHp
+ mHe0 = 4.0 * nHe0
+ mHep = 4.0 * nHep
+ mHepp = 4.0 * nHepp
+ # neglect electron mass contributions
+
+ mtot = mH0 + mHp + mHe0 + mHep + mHepp # + me
+
+ XH0 = mH0 / mtot
+ XHp = mHp / mtot
+ XHe0 = mHe0 / mtot
+ XHep = mHep / mtot
+ XHepp = mHepp / mtot
+ # Xe = me / mtot
+
+ return XH0, XHp, XHe0, XHep, XHepp
+
+
+def internal_energy(T, mu):
+ """
+ Compute the internal energy of the gas for a given
+ temperature and mean molecular weight
+ """
+ # Using u = 1 / (gamma - 1) * p / rho
+ # and p = N/V * kT = rho / (mu * m_u) * kT
+
+ u = unyt.boltzmann_constant * T / (gamma - 1) / (mu * unyt.atomic_mass_unit)
+ return u
+
+
+def mean_molecular_weight(XH0, XHp, XHe0, XHep, XHepp):
+ """
+ Determines the mean molecular weight for given
+ mass fractions of
+ hydrogen: XH0
+ H+: XHp
+ He: XHe0
+ He+: XHep
+ He++: XHepp
+
+ returns:
+ mu: mean molecular weight [in atomic mass units]
+ NOTE: to get the actual mean mass, you still need
+ to multiply it by m_u, as is tradition in the formulae
+ """
+
+ # 1/mu = sum_j X_j / A_j * (1 + E_j)
+ # A_H = 1, E_H = 0
+ # A_Hp = 1, E_Hp = 1
+ # A_He = 4, E_He = 0
+ # A_Hep = 4, E_Hep = 1
+ # A_Hepp = 4, E_Hepp = 2
+ one_over_mu = XH0 + 2 * XHp + 0.25 * XHe0 + 0.5 * XHep + 0.75 * XHepp
+
+ return 1.0 / one_over_mu
+
+
+if __name__ == "__main__":
+
+ glass = h5py.File("glassCube_64.hdf5", "r")
+ parts = glass["PartType0"]
+ xp = parts["Coordinates"][:]
+ h = parts["SmoothingLength"][:]
+ glass.close()
+
+ # find particles closest to the center
+ # and select a couple of them to put the star in their middle
+ r = np.sqrt(np.sum((0.5 - xp) ** 2, axis=1))
+ mininds = np.argsort(r)
+ center_parts = xp[mininds[:4]]
+ xs = center_parts.sum(axis=0) / center_parts.shape[0]
+
+ # Double-check all particles for boundaries
+ for i in range(3):
+ mask = xp[:, i] < 0.0
+ xp[mask, i] += 1.0
+ mask = xp[:, i] > 1.0
+ xp[mask, i] -= 1.0
+
+ # Set up metadata
+ unitL = unyt.Mpc
+ edgelen = 22 * 1e-3 * unitL # 22 so we can cut off 1kpc on each edge for image
+ edgelen = edgelen.to(unitL)
+ boxsize = np.array([1.0, 1.0, 1.0]) * edgelen
+
+ xs = unyt.unyt_array(
+ [np.array([xs[0] * edgelen, xs[1] * edgelen, xs[2] * edgelen])], unitL
+ )
+ xp *= edgelen
+ h *= edgelen
+
+ w = Writer(unit_system=cosmo_units, box_size=boxsize, dimension=3)
+
+ # write particle positions and smoothing lengths
+ w.gas.coordinates = xp
+ w.stars.coordinates = xs
+ w.gas.velocities = np.zeros(xp.shape) * (unitL / unyt.Myr)
+ w.stars.velocities = np.zeros(xs.shape) * (unitL / unyt.Myr)
+ w.gas.smoothing_length = h
+ w.stars.smoothing_length = w.gas.smoothing_length[:1]
+
+ # get gas masses
+ nH = 1e-3 * unyt.cm ** (-3)
+ rhoH = nH * unyt.proton_mass
+ Mtot = rhoH * edgelen ** 3
+ mpart = Mtot / xp.shape[0]
+ mpart = mpart.to(cosmo_units["mass"])
+ w.gas.masses = np.ones(xp.shape[0], dtype=np.float64) * mpart
+ w.stars.masses = np.ones(xs.shape[0], dtype=np.float64) * mpart
+
+ # get gas internal energy for a given temperature and composition
+ XH = 1.0 # hydrogen mass fraction
+ XHe = 0.0 # hydrogen mass fraction
+ T = 100 * unyt.K
+ XHI, XHII, XHeI, XHeII, XHeIII = get_mass_fractions(T, XH, XHe)
+ mu = mean_molecular_weight(XHI, XHII, XHeI, XHeII, XHeIII)
+ internal_energy = internal_energy(T, mu)
+
+ w.gas.internal_energy = np.ones(xp.shape[0], dtype=np.float64) * internal_energy
+
+ w.write("stromgrenSphere-3D.hdf5")
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_3D/makePropagationTestIC.py b/examples/RadiativeTransferTests/StromgrenSphere_3D/makePropagationTestIC.py
new file mode 100755
index 0000000000000000000000000000000000000000..fa81877a54fecb06d24ec848e5bba3e74e4324ec
--- /dev/null
+++ b/examples/RadiativeTransferTests/StromgrenSphere_3D/makePropagationTestIC.py
@@ -0,0 +1,57 @@
+#!/usr/bin/env python3
+
+# ---------------------------------------------------------------------
+# Add a single star in the center of a glass distribution
+# Intended for the photon propagation test.
+# ---------------------------------------------------------------------
+
+from swiftsimio import Writer
+import unyt
+import numpy as np
+import h5py
+
+
+glass = h5py.File("glassCube_64.hdf5", "r")
+parts = glass["PartType0"]
+xp = parts["Coordinates"][:]
+h = parts["SmoothingLength"][:]
+glass.close()
+
+# replace the particle closest to the center
+# by the star
+r = np.sqrt(np.sum((0.5 - xp) ** 2, axis=1))
+rmin = np.argmin(r)
+xs = xp[rmin]
+xp = np.delete(xp, rmin, axis=0)
+h = np.delete(h, rmin)
+
+
+unitL = unyt.cm
+t_end = 1e-3 * unyt.s
+edgelen = unyt.c.to("cm/s") * t_end * 2.0
+edgelen = edgelen.to(unitL)
+boxsize = unyt.unyt_array([edgelen.v, edgelen.v, edgelen.v], unitL)
+
+xs = unyt.unyt_array(
+ [np.array([xs[0] * edgelen, xs[1] * edgelen, xs[2] * edgelen])], unitL
+)
+xp *= edgelen
+h *= edgelen
+
+
+w = Writer(unit_system=unyt.unit_systems.cgs_unit_system, box_size=boxsize, dimension=3)
+
+w.gas.coordinates = xp
+w.stars.coordinates = xs
+w.gas.velocities = np.zeros(xp.shape) * (unyt.cm / unyt.s)
+w.stars.velocities = np.zeros(xs.shape) * (unyt.cm / unyt.s)
+w.gas.masses = np.ones(xp.shape[0], dtype=float) * 100 * unyt.g
+w.stars.masses = np.ones(xs.shape[0], dtype=float) * 100 * unyt.g
+w.gas.internal_energy = (
+ np.ones(xp.shape[0], dtype=float) * (300.0 * unyt.kb * unyt.K) / unyt.g
+)
+
+w.gas.smoothing_length = h
+w.stars.smoothing_length = w.gas.smoothing_length[:1]
+
+w.write("propagationTest-3D.hdf5")
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_3D/plotPhotonPropagationCheck.py b/examples/RadiativeTransferTests/StromgrenSphere_3D/plotPhotonPropagationCheck.py
index 253389a6c8dadc762e3398ea3cde60e903ef8cb0..11516371aa4535c10ff62d82d9f3cad1504fd12f 100755
--- a/examples/RadiativeTransferTests/StromgrenSphere_3D/plotPhotonPropagationCheck.py
+++ b/examples/RadiativeTransferTests/StromgrenSphere_3D/plotPhotonPropagationCheck.py
@@ -16,21 +16,27 @@
# correctly.
# ----------------------------------------------------------------------
-import sys, os, gc
-import swiftsimio
+import gc
+import os
+import sys
+import matplotlib as mpl
import numpy as np
+import swiftsimio
import unyt
from matplotlib import pyplot as plt
-import matplotlib as mpl
from scipy import stats
from scipy.optimize import curve_fit
# Parameters users should/may tweak
# snapshot basename
-snapshot_base = "output"
+snapshot_base = "propagation_test"
+
+# additional anisotropy estimate plot?
+plot_anisotropy_estimate = False
# which photon group to use.
+# NOTE: array index, not group number (which starts at 1 for GEAR)
group_index = 0
scatterplot_kwargs = {
@@ -168,6 +174,7 @@ def plot_photons(filename, emin, emax, fmin, fmax):
time = meta.time
r_expect = meta.time * meta.reduced_lightspeed
+ L = None
use_const_emission_rates = bool(meta.parameters["GEARRT:use_const_emission_rates"])
if use_const_emission_rates:
# read emission rate parameter as string
@@ -186,7 +193,11 @@ def plot_photons(filename, emin, emax, fmin, fmax):
const_emission_rates = unyt.unyt_array(emlist, unyt.L_Sun)
L = const_emission_rates[group_index]
- fig = plt.figure(figsize=(5 * 4, 5.5), dpi=200)
+ if plot_anisotropy_estimate:
+ ncols = 4
+ else:
+ ncols = 3
+ fig = plt.figure(figsize=(5 * ncols, 5.5), dpi=200)
nbins = 100
r_bin_edges = np.linspace(0.5 * edgelen * 1e-3, 0.507 * edgelen, nbins + 1)
@@ -205,8 +216,6 @@ def plot_photons(filename, emin, emax, fmin, fmax):
Fz = getattr(data.gas.photon_fluxes, "Group" + str(group_index + 1) + "Z")
fmag = np.sqrt(Fx ** 2 + Fy ** 2 + Fz ** 2)
- sum_fmag = fmag.sum()
- max_fmag = fmag.max()
particle_count, _ = np.histogram(
r,
bins=r_analytical_bin_edges,
@@ -221,7 +230,7 @@ def plot_photons(filename, emin, emax, fmin, fmax):
# ------------------------
# Plot photon energies
# ------------------------
- ax1 = fig.add_subplot(1, 4, 1)
+ ax1 = fig.add_subplot(1, ncols, 1)
ax1.set_title("Particle Radiation Energies")
ax1.set_ylabel("Photon Energy [$" + energy_units_str + "$]")
@@ -278,7 +287,7 @@ def plot_photons(filename, emin, emax, fmin, fmax):
# ------------------------------
# Plot binned photon energies
# ------------------------------
- ax2 = fig.add_subplot(1, 4, 2)
+ ax2 = fig.add_subplot(1, ncols, 2)
ax2.set_title("Total Radiation Energy in radial bins")
ax2.set_ylabel("Total Photon Energy [$" + energy_units_str + "$]")
@@ -323,7 +332,7 @@ def plot_photons(filename, emin, emax, fmin, fmax):
# ------------------------------
# Plot photon fluxes
# ------------------------------
- ax3 = fig.add_subplot(1, 4, 3)
+ ax3 = fig.add_subplot(1, ncols, 3)
ax3.set_title("Particle Radiation Flux Magnitudes")
ax3.set_ylabel("Photon Flux Magnitude [$" + flux_units_str + "$]")
@@ -387,55 +396,58 @@ def plot_photons(filename, emin, emax, fmin, fmax):
# Plot photon flux sum
# ------------------------------
- ax4 = fig.add_subplot(1, 4, 4)
- ax4.set_title("Vectorial Sum of Radiation Flux in radial bins")
- ax4.set_ylabel("[1]")
-
- fmag_sum_bin, _, _ = stats.binned_statistic(
- r,
- fmag,
- statistic="sum",
- bins=r_bin_edges,
- range=(r_bin_edges[0], r_bin_edges[-1]),
- )
- mask_sum = fmag_sum_bin > 0
- fmag_max_bin, _, _ = stats.binned_statistic(
- r,
- fmag,
- statistic="max",
- bins=r_bin_edges,
- range=(r_bin_edges[0], r_bin_edges[-1]),
- )
- mask_max = fmag_max_bin > 0
- Fx_sum_bin, _, _ = stats.binned_statistic(
- r,
- Fx,
- statistic="sum",
- bins=r_bin_edges,
- range=(r_bin_edges[0], r_bin_edges[-1]),
- )
- Fy_sum_bin, _, _ = stats.binned_statistic(
- r,
- Fy,
- statistic="sum",
- bins=r_bin_edges,
- range=(r_bin_edges[0], r_bin_edges[-1]),
- )
- F_sum_bin = np.sqrt(Fx_sum_bin ** 2 + Fy_sum_bin ** 2)
+ if plot_anisotropy_estimate:
+ ax4 = fig.add_subplot(1, ncols, 4)
+ ax4.set_title("Vectorial Sum of Radiation Flux in radial bins")
+ ax4.set_ylabel("[1]")
+
+ fmag_sum_bin, _, _ = stats.binned_statistic(
+ r,
+ fmag,
+ statistic="sum",
+ bins=r_bin_edges,
+ range=(r_bin_edges[0], r_bin_edges[-1]),
+ )
+ mask_sum = fmag_sum_bin > 0
+ fmag_max_bin, _, _ = stats.binned_statistic(
+ r,
+ fmag,
+ statistic="max",
+ bins=r_bin_edges,
+ range=(r_bin_edges[0], r_bin_edges[-1]),
+ )
+ mask_max = fmag_max_bin > 0
+ Fx_sum_bin, _, _ = stats.binned_statistic(
+ r,
+ Fx,
+ statistic="sum",
+ bins=r_bin_edges,
+ range=(r_bin_edges[0], r_bin_edges[-1]),
+ )
+ Fy_sum_bin, _, _ = stats.binned_statistic(
+ r,
+ Fy,
+ statistic="sum",
+ bins=r_bin_edges,
+ range=(r_bin_edges[0], r_bin_edges[-1]),
+ )
+ F_sum_bin = np.sqrt(Fx_sum_bin ** 2 + Fy_sum_bin ** 2)
- ax4.plot(
- r_bin_centres[mask_sum],
- F_sum_bin[mask_sum] / fmag_sum_bin[mask_sum],
- **lineplot_kwargs,
- label="$\left| \sum_{i \in \mathrm{particles \ in \ bin}} \mathbf{F}_i \\right| $ / $\sum_{i \in \mathrm{particles \ in \ bin}} \left| \mathbf{F}_{i} \\right| $",
- )
- ax4.plot(
- r_bin_centres[mask_max],
- F_sum_bin[mask_max] / fmag_max_bin[mask_max],
- **lineplot_kwargs,
- linestyle="--",
- label="$\left| \sum_{i \in \mathrm{particles \ in \ bin}} \mathbf{F}_i \\right| $ / $\max_{i \in \mathrm{particles \ in \ bin}} \left| \mathbf{F}_{i} \\right| $",
- )
+ ax4.plot(
+ r_bin_centres[mask_sum],
+ F_sum_bin[mask_sum] / fmag_sum_bin[mask_sum],
+ **lineplot_kwargs,
+ label="$\left| \sum_{i \in \mathrm{particles \ in \ bin}} \mathbf{F}_i \\right| $ "
+ + "/ $\sum_{i \in \mathrm{particles \ in \ bin}} \left| \mathbf{F}_{i} \\right| $",
+ )
+ ax4.plot(
+ r_bin_centres[mask_max],
+ F_sum_bin[mask_max] / fmag_max_bin[mask_max],
+ **lineplot_kwargs,
+ linestyle="--",
+ label="$\left| \sum_{i \in \mathrm{particles \ in \ bin}} \mathbf{F}_i \\right| $ "
+ + " / $\max_{i \in \mathrm{particles \ in \ bin}} \left| \mathbf{F}_{i} \\right| $",
+ )
# -------------------------------------------
# Cosmetics that all axes have in common
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_3D/plotRadiationProjection.py b/examples/RadiativeTransferTests/StromgrenSphere_3D/plotRadiationProjection.py
index 05cfa4a6a82dda75934ebeb2914c95a2129d0a6c..617226cc71433109cba19160769b7fc219882fd2 100755
--- a/examples/RadiativeTransferTests/StromgrenSphere_3D/plotRadiationProjection.py
+++ b/examples/RadiativeTransferTests/StromgrenSphere_3D/plotRadiationProjection.py
@@ -23,8 +23,10 @@ from matplotlib.colors import SymLogNorm
# plot all groups and all photon quantities
plot_all_data = True
-# snapshot basename
-snapshot_base = "output"
+
+# plotting propagation tests or Stromgren sphere?
+do_stromgren_sphere = True
+
# fancy up the plots a bit?
fancy = True
@@ -34,12 +36,24 @@ imshow_kwargs = {"origin": "lower", "cmap": "viridis"}
# parameters for swiftsimio projections
projection_kwargs = {"resolution": 1024, "parallel": True}
+# snapshot basename
+snapshot_base = "propagation_test"
+
# Set Units of your choice
energy_units = unyt.erg
energy_units_str = "\\rm{erg}"
flux_units = 1e10 * energy_units / unyt.cm ** 2 / unyt.s
flux_units_str = "10^{10} \\rm{erg} \\ \\rm{cm}^{-2} \\ \\rm{s}^{-1}"
time_units = unyt.s
+
+if do_stromgren_sphere:
+ snapshot_base = "output"
+
+ energy_units = 1e50 * unyt.erg
+ energy_units_str = "10^{50} \\rm{erg}"
+ flux_units = 1e50 * unyt.erg / unyt.kpc ** 2 / unyt.Gyr
+ flux_units_str = "10^{60} \\rm{erg} \\ \\rm{kpc}^{-2} \\ \\rm{Gyr}^{-1}"
+ time_units = unyt.Myr
# -----------------------------------------------------------------------
@@ -76,6 +90,10 @@ def get_snapshot_list(snapshot_basename="output"):
quit(1)
snaplist.append(fname)
+ if len(snaplist) == 0:
+ print("Didn't find any snapshots with basename '", snapshot_basename, "'")
+ quit(1)
+
return snaplist
@@ -142,7 +160,7 @@ def plot_photons(filename, energy_boundaries=None, flux_boundaries=None):
if plot_all_data:
fig = plt.figure(figsize=(5 * 4, 5.05 * ngroups), dpi=200)
- figname = filename[:-5] + "-all-quantities.png"
+ figname = filename[:-5] + "-radiation-projection.png"
for g in range(ngroups):
@@ -320,7 +338,6 @@ def get_minmax_vals(snaplist):
dirmin = []
dirmax = []
for direction in ["X", "Y", "Z"]:
- new_attribute_str = "radiation_flux" + str(g + 1) + direction
f = getattr(data.gas.photon_fluxes, "Group" + str(g + 1) + direction)
dirmin.append(f.min())
dirmax.append(f.max())
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_3D/plotSolution.py b/examples/RadiativeTransferTests/StromgrenSphere_3D/plotSolution.py
new file mode 100755
index 0000000000000000000000000000000000000000..1c00fb71c14f94583e8a9b4532bd91f5339bc474
--- /dev/null
+++ b/examples/RadiativeTransferTests/StromgrenSphere_3D/plotSolution.py
@@ -0,0 +1,285 @@
+#!/usr/bin/env python3
+
+# ----------------------------------------------------
+# Plot slices of the hydrogen number density,
+# pressure, temperature, and hydrogen mass fraction.
+# ----------------------------------------------------
+
+import sys
+import os
+import swiftsimio
+import gc
+import unyt
+from matplotlib import pyplot as plt
+import matplotlib as mpl
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+from matplotlib.colors import LogNorm
+from swiftsimio.visualisation.slice import slice_gas
+
+# Parameters users should/may tweak
+
+# snapshot basename
+snapshot_base = "output"
+
+# parameters for imshow plots
+imshow_kwargs = {"origin": "lower"}
+
+# parameters for swiftsimio slices
+slice_kwargs = {"slice": 0.5, "resolution": 1000, "parallel": 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 get_snapshot_list(snapshot_basename="output"):
+ """
+ Find the snapshot(s) that are to be plotted
+ and return their names as list
+ """
+
+ snaplist = []
+
+ if plot_all:
+ dirlist = os.listdir()
+ for f in dirlist:
+ if f.startswith(snapshot_basename) and f.endswith("hdf5"):
+ snaplist.append(f)
+
+ snaplist = sorted(snaplist)
+
+ else:
+ fname = snapshot_basename + "_" + str(snapnr).zfill(4) + ".hdf5"
+ if not os.path.exists(fname):
+ print("Didn't find file", fname)
+ quit(1)
+ snaplist.append(fname)
+
+ if len(snaplist) == 0:
+ print("Didn't find any outputs with basename '", snapshot_basename, "'")
+ quit(1)
+
+ return snaplist
+
+
+def mean_molecular_weight(XH0, XHp, XHe0, XHep, XHepp):
+ """
+ Determines the mean molecular weight for given
+ mass fractions of
+ hydrogen: XH0
+ H+: XHp
+ He: XHe0
+ He+: XHep
+ He++: XHepp
+
+ returns:
+ mu: mean molecular weight [in atomic mass units]
+ NOTE: to get the actual mean mass, you still need
+ to multiply it by m_u, as is tradition in the formulae
+ """
+
+ # 1/mu = sum_j X_j / A_j * (1 + E_j)
+ # A_H = 1, E_H = 0
+ # A_Hp = 1, E_Hp = 1
+ # A_He = 4, E_He = 0
+ # A_Hep = 4, E_Hep = 1
+ # A_Hepp = 4, E_Hepp = 2
+ one_over_mu = XH0 + 2 * XHp + 0.25 * XHe0 + 0.5 * XHep + 0.75 * XHepp
+
+ return 1.0 / one_over_mu
+
+
+def gas_temperature(u, mu, gamma):
+ """
+ Compute the gas temperature given the specific internal
+ energy u and the mean molecular weight mu
+ """
+
+ # Using u = 1 / (gamma - 1) * p / rho
+ # and p = N/V * kT = rho / (mu * m_u) * kT
+
+ T = u * (gamma - 1) * mu * unyt.atomic_mass_unit / unyt.boltzmann_constant
+
+ return T.to("K")
+
+
+def set_colorbar(ax, im):
+ """
+ Adapt the colorbar a bit for axis object <ax> and
+ imshow instance <im>
+ """
+ divider = make_axes_locatable(ax)
+ cax = divider.append_axes("right", size="5%", pad=0.05)
+ plt.colorbar(im, cax=cax)
+ return
+
+
+def plot_result(filename):
+ """
+ Create and save the plot
+ """
+ print("working on", filename)
+
+ data = swiftsimio.load(filename)
+ meta = data.metadata
+
+ global imshow_kwargs
+ imshow_kwargs["extent"] = [
+ 0.0 * meta.boxsize[0].v,
+ 0.9 * meta.boxsize[0].v,
+ 0.0 * meta.boxsize[1].v,
+ 0.9 * meta.boxsize[1].v,
+ ]
+ cutoff = int(0.05 * slice_kwargs["resolution"])
+
+ mass_map = slice_gas(data, project="masses", **slice_kwargs)
+ gamma = meta.hydro_scheme["Adiabatic index"][0]
+
+ data.gas.mXHI = data.gas.ion_mass_fractions.HI * data.gas.masses
+ data.gas.mXHII = data.gas.ion_mass_fractions.HII * data.gas.masses
+ data.gas.mP = data.gas.pressures * data.gas.masses
+ data.gas.mrho = data.gas.densities * data.gas.masses
+
+ imf = data.gas.ion_mass_fractions
+ mu = mean_molecular_weight(imf.HI, imf.HII, imf.HeI, imf.HeII, imf.HeIII)
+ data.gas.mT = (
+ gas_temperature(data.gas.internal_energies, mu, gamma) * data.gas.masses
+ )
+
+ mass_weighted_hydrogen_map = slice_gas(data, project="mXHI", **slice_kwargs)
+ mass_weighted_pressure_map = slice_gas(data, project="mP", **slice_kwargs)
+ mass_weighted_density_map = slice_gas(data, project="mrho", **slice_kwargs)
+ mass_weighted_temperature_map = slice_gas(data, project="mT", **slice_kwargs)
+
+ hydrogen_map = mass_weighted_hydrogen_map / mass_map
+ hydrogen_map = hydrogen_map[cutoff:-cutoff, cutoff:-cutoff]
+
+ pressure_map = mass_weighted_pressure_map / mass_map
+ pressure_map = pressure_map[cutoff:-cutoff, cutoff:-cutoff]
+ pressure_map = pressure_map.to("g/cm/s**2")
+
+ density_map = mass_weighted_density_map / mass_map
+ density_map = density_map[cutoff:-cutoff, cutoff:-cutoff]
+ density_map = density_map.to("kg/cm**3")
+ density_map = density_map / unyt.proton_mass
+
+ temperature_map = mass_weighted_temperature_map / mass_map
+ temperature_map = temperature_map[cutoff:-cutoff, cutoff:-cutoff]
+ temperature_map = temperature_map.to("K")
+
+ fig = plt.figure(figsize=(12, 12), dpi=200)
+ figname = filename[:-5] + ".png"
+
+ ax1 = fig.add_subplot(221)
+ ax2 = fig.add_subplot(222)
+ ax3 = fig.add_subplot(223)
+ ax4 = fig.add_subplot(224)
+
+ try:
+ im1 = ax1.imshow(
+ density_map.T,
+ **imshow_kwargs,
+ norm=LogNorm(vmin=1e-4, vmax=1e-1),
+ cmap="bone",
+ )
+ set_colorbar(ax1, im1)
+ ax1.set_title(r"Hydrogen Number Density [cm$^{-3}$]")
+ except ValueError:
+ print(
+ filename,
+ "densities wrong? min",
+ data.gas.densities.min(),
+ "max",
+ data.gas.densities.max(),
+ )
+ return
+
+ try:
+ im2 = ax2.imshow(
+ hydrogen_map.T,
+ **imshow_kwargs,
+ norm=LogNorm(vmin=1e-3, vmax=1.0),
+ cmap="cividis",
+ )
+ set_colorbar(ax2, im2)
+ ax2.set_title("Hydrogen Mass Fraction [1]")
+ except ValueError:
+ print(
+ filename,
+ "mass fraction wrong? min",
+ data.gas.ion_mass_fractions.HI.min(),
+ "max",
+ data.gas.ion_mass_fractions.HI.max(),
+ )
+ return
+
+ try:
+ im3 = ax3.imshow(
+ pressure_map.T,
+ **imshow_kwargs,
+ norm=LogNorm(vmin=1e-15, vmax=1e-12),
+ cmap="viridis",
+ )
+ set_colorbar(ax3, im3)
+ ax3.set_title(r"Pressure [g/cm/s$^2$]")
+ except ValueError:
+ print(
+ filename,
+ "pressures wrong? min",
+ data.gas.pressures.min(),
+ "max",
+ data.gas.pressures.max(),
+ )
+ return
+
+ try:
+ im4 = ax4.imshow(
+ temperature_map.T,
+ **imshow_kwargs,
+ norm=LogNorm(vmin=1e2, vmax=4e4),
+ cmap="inferno",
+ )
+ set_colorbar(ax4, im4)
+ ax4.set_title(r"Temperature [K]")
+ except ValueError:
+ print(
+ filename,
+ "temperatures wrong? min",
+ temperature_map.min(),
+ "max",
+ temperature_map.max(),
+ )
+ return
+
+ for ax in [ax1, ax2, ax3, ax4]:
+ ax.set_xlabel("[kpc]")
+ ax.set_ylabel("[kpc]")
+
+ title = filename.replace("_", "\_") # exception handle underscore for latex
+ if meta.cosmology is not None:
+ title += ", $z$ = {0:.2e}".format(meta.z)
+ title += ", $t$ = {0:.2e}".format(meta.time.to("Myr"))
+ fig.suptitle(title)
+
+ plt.tight_layout()
+ plt.savefig(figname)
+ plt.close()
+ gc.collect()
+ return
+
+
+if __name__ == "__main__":
+
+ snaplist = get_snapshot_list(snapshot_base)
+ # snaplist = snaplist[118:]
+
+ for f in snaplist:
+ plot_result(f)
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_3D/propagationTest-3D.yml b/examples/RadiativeTransferTests/StromgrenSphere_3D/propagationTest-3D.yml
new file mode 100644
index 0000000000000000000000000000000000000000..b47058bbb49e4c396b83c4794ff1fc9c562ca3e1
--- /dev/null
+++ b/examples/RadiativeTransferTests/StromgrenSphere_3D/propagationTest-3D.yml
@@ -0,0 +1,60 @@
+MetaData:
+ run_name: StromgrenSpherePropagationTest3D
+
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.
+ UnitLength_in_cgs: 1.
+ UnitVelocity_in_cgs: 1.
+ UnitCurrent_in_cgs: 1.
+ UnitTemp_in_cgs: 1.
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 0.001 # 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: 2.e-04 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: propagation_test # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 0.0001 # Time between snapshots
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ time_first: 0.
+ delta_time: 5e-4 # 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.6 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 10. # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./propagationTest-3D.hdf5 # The file to read
+ periodic: 1 # peridioc ICs. Keep them periodic so we don't loose photon energy. TODO: CHANGE LATER WHEN YOU ACTUALLY DO GAS INTERACTIONS
+
+Scheduler:
+ max_top_level_cells: 64
+
+GEARRT:
+ f_reduce_c: 1. # reduce the speed of light for the RT solver by multiplying c with this factor
+ CFL_condition: 0.99
+ photon_groups_Hz: [] # Photon frequency group bin edges in Hz. Needs to be 1 less than the number of groups (N) requested during the configuration (--with-RT=GEAR_N).
+ use_const_emission_rates: 1 # (Optional) use constant emission rates for stars as defined with star_emission_rates_erg_LSol parameter
+ star_emission_rates_LSol: [1e-28] # (Optional) constant star emission rates for each photon frequency group to use if use_constant_emission_rates is set.
+ hydrogen_mass_fraction: 0.76 # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
+ stellar_spectrum_type: 0 # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
+ stellar_spectrum_const_max_frequency_Hz: 1.e17 # (Conditional) if stellar_spectrum_type=0, use this maximal frequency for the constant spectrum.
+ set_initial_ionization_mass_fractions: 1 # (Optional) manually overwrite initial mass fraction of each species (using the values you set below)
+ mass_fraction_HI: 0.76 # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HI mass fractions to this value
+ mass_fraction_HII: 0. # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HII mass fractions to this value
+ mass_fraction_HeI: 0.24 # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeI mass fractions to this value
+ mass_fraction_HeII: 0. # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeII mass fractions to this value
+ mass_fraction_HeIII: 0. # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeIII mass fractions to this value
+ skip_thermochemistry: 1 # skip thermochemistry.
+ stars_max_timestep: 1.562500e-05 # (Optional) restrict the maximal timestep of stars to this value (in internal units). Set to negative to turn off.
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_3D/run.sh b/examples/RadiativeTransferTests/StromgrenSphere_3D/run.sh
index 0afaf94220f0aa102c5db22c05c79e47ec2da6e0..1e60f3e825a4cecfdd7540d371421e0a6970d05b 100755
--- a/examples/RadiativeTransferTests/StromgrenSphere_3D/run.sh
+++ b/examples/RadiativeTransferTests/StromgrenSphere_3D/run.sh
@@ -24,5 +24,4 @@ fi
# Plot the photon propagation checks.
# Make sure you set the correct photon group to plot
# inside the script
-python3 ./plotPhotonPropagationCheck.py
-
+python3 ./plotSolution.py
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_3D/runPropagationTest.sh b/examples/RadiativeTransferTests/StromgrenSphere_3D/runPropagationTest.sh
new file mode 100755
index 0000000000000000000000000000000000000000..c61fc9de7dd0badf1efcf26f0f73ecf9a872b521
--- /dev/null
+++ b/examples/RadiativeTransferTests/StromgrenSphere_3D/runPropagationTest.sh
@@ -0,0 +1,27 @@
+#!/bin/bash
+
+# make run.sh fail if a subcommand fails
+set -e
+set -o pipefail
+
+if [ ! -e glassCube_64.hdf5 ]
+then
+ echo "Fetching initial glass file for Strömgen Sphere 3D example ..."
+ ./getGlass.sh
+fi
+
+if [ ! -f 'propagationTest-3D.hdf5' ]; then
+ echo "Generating ICs"
+ python3 makePropagationTestIC.py
+fi
+
+# Run SWIFT with RT
+../../swift \
+ --hydro --threads=4 --stars --external-gravity \
+ --feedback --radiation --fpe \
+ ./propagationTest-3D.yml 2>&1 | tee output.log
+
+# Plot the photon propagation checks.
+# Make sure you set the correct photon group to plot
+# inside the script
+python3 ./plotPhotonPropagationCheck.py
diff --git a/examples/RadiativeTransferTests/StromgrenSphere_3D/stromgrenSphere-3D.yml b/examples/RadiativeTransferTests/StromgrenSphere_3D/stromgrenSphere-3D.yml
index 97499f58f4342e14c827b97d062a9129247d1484..048c948d61a88dfb4cd5e9bce23db7dc524d2ff6 100644
--- a/examples/RadiativeTransferTests/StromgrenSphere_3D/stromgrenSphere-3D.yml
+++ b/examples/RadiativeTransferTests/StromgrenSphere_3D/stromgrenSphere-3D.yml
@@ -3,57 +3,62 @@ MetaData:
# Define the system of units to use internally.
InternalUnitSystem:
- UnitMass_in_cgs: 1.
- UnitLength_in_cgs: 1.
- UnitVelocity_in_cgs: 1.
+ 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.
+ UnitTemp_in_cgs: 1. # K
+
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 0.001 # end time: radiation reaches edge of box
+ time_end: 0.512 # 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: 2.e-04 # 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).
# Parameters governing the snapshots
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.0001 # Time between snapshots
+ delta_time: 0.001 # Time between snapshots
# Parameters governing the conserved quantities statistics
Statistics:
time_first: 0.
- delta_time: 5e-4 # Time between statistics output
+ 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.6 # Courant-Friedrich-Levy condition for time integration.
+ CFL_condition: 0.8 # Courant-Friedrich-Levy condition for time integration.
minimal_temperature: 10. # Kelvin
# Parameters related to the initial conditions
InitialConditions:
file_name: ./stromgrenSphere-3D.hdf5 # The file to read
- periodic: 1 # peridioc ICs. Keep them periodic so we don't loose photon energy. TODO: CHANGE LATER WHEN YOU ACTUALLY DO GAS INTERACTIONS
+ periodic: 0 # periodic ICs.
Scheduler:
max_top_level_cells: 64
+Stars:
+ resolution_eta: 2.2348 # (Optional) Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). Defaults to the SPH value.
+
GEARRT:
- f_reduce_c: 1. # reduce the speed of light for the RT solver by multiplying c with this factor
- CFL_condition: 0.6
- photon_groups_Hz: [] # Photon frequency group bin edges in Hz. Needs to be 1 less than the number of groups (N) requested during the configuration (--with-RT=GEAR_N).
- use_const_emission_rates: 1 # (Optional) use constant emission rates for stars as defined with star_emission_rates_erg_LSol parameter
- star_emission_rates_LSol: [1e-28] # (Optional) constant star emission rates for each photon frequency group to use if use_constant_emission_rates is set.
- hydrogen_mass_fraction: 0.76 # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
- stellar_spectrum_type: 0 # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
- stellar_spectrum_const_max_frequency_Hz: 1.e17 # (Conditional) if stellar_spectrum_type=0, use this maximal frequency for the constant spectrum.
+ f_reduce_c: 0.01 # reduce the speed of light for the RT solver by multiplying c with this factor
+ CFL_condition: 0.9
+ photon_groups_Hz: [3.288e15, 5.945e15, 13.157e15] # Photon frequency group bin edges in Hz. Needs to be 1 less than the number of groups (N) requested during the configuration (--with-RT=GEAR_N). Outer edges of zero and infinity are assumed.
+ use_const_emission_rates: 1 # (Optional) use constant emission rates for stars as defined with star_emission_rates_erg_LSol parameter
+ star_emission_rates_LSol: [247985.76636, 247985.76636] # (Optional) constant star emission rates for each photon frequency group to use if use_constant_emission_rates is set, in units of Solar Luminosity.
+ hydrogen_mass_fraction: 1.00 # 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: 1 # (Optional) manually overwrite initial mass fraction of each species (using the values you set below)
- mass_fraction_HI: 0.76 # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HI mass fractions to this value
- mass_fraction_HII: 0. # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HII mass fractions to this value
- mass_fraction_HeI: 0.24 # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeI mass fractions to this value
+ mass_fraction_HI: 0.999999 # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HI mass fractions to this value
+ mass_fraction_HII: 1.e-6 # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HII mass fractions to this value
+ mass_fraction_HeI: 0. # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeI mass fractions to this value
mass_fraction_HeII: 0. # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeII mass fractions to this value
mass_fraction_HeIII: 0. # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeIII mass fractions to this value
-
+ stellar_spectrum_type: 1 # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 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.
+ stars_max_timestep: 5e-4 # (Optional) restrict the maximal timestep of stars to this value (in internal units)
diff --git a/examples/RadiativeTransferTests/UniformBox_3D/README b/examples/RadiativeTransferTests/UniformBox_3D/README
index 6305fd99428c497843e590d218b6b6326632999b..dd27a9062be5d36c8a615a9d29b2ba44e2b07ff7 100644
--- a/examples/RadiativeTransferTests/UniformBox_3D/README
+++ b/examples/RadiativeTransferTests/UniformBox_3D/README
@@ -29,4 +29,3 @@ Tl;dr:
GEAR scheme with debugging checks enabled
- `./rt_sanity_checks-GEAR.py` is made to work on any run with the GEAR
scheme and debugging checks enabled.
-
diff --git a/examples/RadiativeTransferTests/UniformBox_3D/makeIC.py b/examples/RadiativeTransferTests/UniformBox_3D/makeIC.py
index bcfd67f510f2a89cf15940f8b68cac5fc8aad7b5..ff47c49a102869c399d55b2d62455be63839c288 100755
--- a/examples/RadiativeTransferTests/UniformBox_3D/makeIC.py
+++ b/examples/RadiativeTransferTests/UniformBox_3D/makeIC.py
@@ -53,7 +53,7 @@ for i in range(n_p):
y = (j + 0.5) * dx
for k in range(n_p):
z = (k + 0.5) * dx
- xp.append(np.array([x, y, z], dtype=np.float))
+ xp.append(np.array([x, y, z], dtype=np.float64))
# Generate star coordinates
for i in range(n_s):
@@ -63,7 +63,7 @@ for i in range(n_s):
y = 0.4 * boxsize + (j + 0.52) * ds
for k in range(n_s):
z = 0.4 * boxsize + (k + 0.52) * ds
- xs.append(np.array([x, y, z], dtype=np.float))
+ xs.append(np.array([x, y, z], dtype=np.float64))
xp = unyt.unyt_array(xp, boxsize.units)
xs = unyt.unyt_array(xs, boxsize.units)
@@ -75,10 +75,10 @@ w.gas.coordinates = xp
w.stars.coordinates = xs
w.gas.velocities = np.zeros(xp.shape) * (unyt.cm / unyt.s)
w.stars.velocities = np.zeros(xs.shape) * (unyt.cm / unyt.s)
-w.gas.masses = np.ones(xp.shape[0], dtype=np.float) * 1000 * unyt.g
-w.stars.masses = np.ones(xs.shape[0], dtype=np.float) * 1000 * unyt.g
+w.gas.masses = np.ones(xp.shape[0], dtype=np.float64) * 1000 * unyt.g
+w.stars.masses = np.ones(xs.shape[0], dtype=np.float64) * 1000 * unyt.g
w.gas.internal_energy = (
- np.ones(xp.shape[0], dtype=np.float) * (300.0 * unyt.kb * unyt.K) / (unyt.g)
+ np.ones(xp.shape[0], dtype=np.float64) * (300.0 * unyt.kb * unyt.K) / unyt.g
)
# Generate initial guess for smoothing lengths based on MIPS
@@ -110,12 +110,11 @@ parts = F["/PartType0"]
for grp in range(nPhotonGroups):
dsetname = "PhotonEnergiesGroup{0:d}".format(grp + 1)
- energydata = np.ones((nparts), dtype=np.float) * (grp + 1)
+ energydata = np.ones((nparts), dtype=np.float64) * (grp + 1)
parts.create_dataset(dsetname, data=energydata)
dsetname = "PhotonFluxesGroup{0:d}".format(grp + 1)
- # if dsetname not in parts.keys():
- fluxdata = np.ones((nparts, 3), dtype=np.float) * (grp + 1)
+ fluxdata = np.ones((nparts, 3), dtype=np.float64) * (grp + 1)
fluxdata[:, 1] *= 2.0
fluxdata[:, 2] *= 3.0
parts.create_dataset(dsetname, data=fluxdata)
diff --git a/examples/RadiativeTransferTests/UniformBox_3D/plotRadiationProjection.py b/examples/RadiativeTransferTests/UniformBox_3D/plotRadiationProjection.py
index aec820ea3a62b8abd314f95d3b7f7ebc6d493581..c820b283139ba8eb527d190b4e07a89fadb9b18b 100755
--- a/examples/RadiativeTransferTests/UniformBox_3D/plotRadiationProjection.py
+++ b/examples/RadiativeTransferTests/UniformBox_3D/plotRadiationProjection.py
@@ -320,7 +320,6 @@ def get_minmax_vals(snaplist):
dirmin = []
dirmax = []
for direction in ["X", "Y", "Z"]:
- new_attribute_str = "radiation_flux" + str(g + 1) + direction
f = getattr(data.gas.photon_fluxes, "Group" + str(g + 1) + direction)
dirmin.append(f.min())
dirmax.append(f.max())
diff --git a/examples/RadiativeTransferTests/UniformBox_3D/rt_sanity_checks-GEAR.py b/examples/RadiativeTransferTests/UniformBox_3D/rt_sanity_checks-GEAR.py
index 62392fe87b478bcda916ff8523e903a5ac300cbe..a9af4f62b94dfe3563fa3145979be1d05e32d895 100755
--- a/examples/RadiativeTransferTests/UniformBox_3D/rt_sanity_checks-GEAR.py
+++ b/examples/RadiativeTransferTests/UniformBox_3D/rt_sanity_checks-GEAR.py
@@ -124,14 +124,12 @@ def check_injection(snapdata, rundata):
ngroups = rundata.ngroups
initial_energies = np.sum(snapdata[0].gas.PhotonEnergies, axis=0)
- initial_time = snapdata[0].time
# Check 1a) : sum initial energy + sum injected energy = sum current energy
# --------------------------------------------------------------------------
# TODO: this assumes no cosmological expansion
for snap in snapdata[1:]:
- dt = snap.time - initial_time
# sum of each group over all particles
photon_energies = np.sum(snap.gas.PhotonEnergies, axis=0)
if snap.has_stars:
@@ -140,7 +138,7 @@ def check_injection(snapdata, rundata):
np.sum(snap.stars.InjectedPhotonEnergy, axis=0)
)
else:
- injected_energies = np.zeros(ngroups, dtype=np.float)
+ injected_energies = np.zeros(ngroups, dtype=np.float64)
for g in range(ngroups):
energy_expected = initial_energies[g] + injected_energies[g]
@@ -183,7 +181,9 @@ def check_injection(snapdata, rundata):
if snapdata[0].has_stars:
emission_at_initial_time = snapdata[0].stars.InjectedPhotonEnergy.sum(axis=0)
else:
- emission_at_initial_time = np.zeros(rundata.ngroups, dtype=np.float) * unyt.erg
+ emission_at_initial_time = (
+ np.zeros(rundata.ngroups, dtype=np.float64) * unyt.erg
+ )
if rundata.use_const_emission_rate and not rundata.hydro_controlled_injection:
if len(snapdata) <= 2:
diff --git a/examples/RadiativeTransferTests/UniformBox_3D/rt_sanity_checks.py b/examples/RadiativeTransferTests/UniformBox_3D/rt_sanity_checks.py
index 6214eba725a34adc1eefe1956b463c1db16fab73..86eafab768de3b1482bc6fdbb6fd4b357a6f6981 100755
--- a/examples/RadiativeTransferTests/UniformBox_3D/rt_sanity_checks.py
+++ b/examples/RadiativeTransferTests/UniformBox_3D/rt_sanity_checks.py
@@ -67,7 +67,6 @@ def check_hydro_sanity(snapdata):
- RT transport calls >= RT gradient calls?
"""
- nsnaps = len(snapdata)
npart = snapdata[0].gas.coords.shape[0]
print("Checking hydro")
@@ -78,7 +77,7 @@ def check_hydro_sanity(snapdata):
for snap in snapdata:
gas = snap.gas
- stars = snap.stars
+ # stars = snap.stars
# has a particle been called at least once?
called = (
@@ -221,8 +220,6 @@ def check_stars_sanity(snapdata):
print("Checking stars")
- nsnaps = len(snapdata)
-
# ----------------------------------------------
# check consistency of individual snapshots
# ----------------------------------------------
@@ -260,9 +257,6 @@ def check_stars_hydro_interaction_sanity(snapdata):
- are total calls each step equal?
"""
- nsnaps = len(snapdata)
- npart = snapdata[0].gas.coords.shape[0]
-
print("Checking hydro vs stars")
# ----------------------------------------------
diff --git a/examples/RadiativeTransferTests/UniformBox_3D/rt_uniform_box_checks-GEAR.py b/examples/RadiativeTransferTests/UniformBox_3D/rt_uniform_box_checks-GEAR.py
index d521a77d79e72ef531055a7914c2cdab89fcef36..0d855cc6529cd59c0a785726e0be039c847592a2 100755
--- a/examples/RadiativeTransferTests/UniformBox_3D/rt_uniform_box_checks-GEAR.py
+++ b/examples/RadiativeTransferTests/UniformBox_3D/rt_uniform_box_checks-GEAR.py
@@ -41,7 +41,6 @@
import numpy as np
-import unyt
from sys import argv
from swift_rt_GEAR_io import get_snap_data
diff --git a/examples/RadiativeTransferTests/UniformBox_3D/swift_rt_GEAR_io.py b/examples/RadiativeTransferTests/UniformBox_3D/swift_rt_GEAR_io.py
index f1e6d88836c07756e8e5fed489c5c5b8affcba9a..623653ba02d0a6bfe5f9974f404a36c85dad4c98 100644
--- a/examples/RadiativeTransferTests/UniformBox_3D/swift_rt_GEAR_io.py
+++ b/examples/RadiativeTransferTests/UniformBox_3D/swift_rt_GEAR_io.py
@@ -26,8 +26,6 @@
import os
-import h5py
-import numpy as np
import unyt
import swiftsimio
@@ -166,7 +164,6 @@ def get_snap_data(prefix="output", skip_snap_zero=False, skip_last_snap=False):
"These tests only work for the GEAR RT scheme.",
"Compile swift --with-rt=GEAR_N",
)
- quit()
ngroups = int(firstfile.metadata.subgrid_scheme["PhotonGroupNumber"])
rundata.ngroups = ngroups
@@ -266,7 +263,7 @@ def get_snap_data(prefix="output", skip_snap_zero=False, skip_last_snap=False):
newsnap.stars = Stars
newsnap.nstars = Stars.IDs.shape[0]
except AttributeError:
- Stars = RTStarData()
+ newsnap.stars = RTStarData()
newsnap.has_stars = False
newsnap.nstars = 0
diff --git a/src/rt/GEAR/rt.h b/src/rt/GEAR/rt.h
index 90350dffaa6108495b894b78ff82b4cd78324b21..d15f601c12ca239aea3aeb5bba8f6310af833ca4 100644
--- a/src/rt/GEAR/rt.h
+++ b/src/rt/GEAR/rt.h
@@ -467,12 +467,13 @@ __attribute__((always_inline)) INLINE static void rt_finalise_transport(
struct rt_part_data* restrict rtd = &p->rt_data;
for (int g = 0; g < RT_NGROUPS; g++) {
+ const float e_old = rtd->conserved[g].energy;
rtd->conserved[g].energy += rtd->flux[g].energy * dt;
rtd->conserved[g].flux[0] += rtd->flux[g].flux[0] * dt;
rtd->conserved[g].flux[1] += rtd->flux[g].flux[1] * dt;
rtd->conserved[g].flux[2] += rtd->flux[g].flux[2] * dt;
rt_check_unphysical_conserved(&rtd->conserved[g].energy,
- rtd->conserved[g].flux, 1);
+ rtd->conserved[g].flux, e_old, 1);
}
}
diff --git a/src/rt/GEAR/rt_getters.h b/src/rt/GEAR/rt_getters.h
index 5c04f4755c6b05812535863938148ee45381f284..c4c2d3f8795c2371c87a1e591181cb03d02e84a8 100644
--- a/src/rt/GEAR/rt_getters.h
+++ b/src/rt/GEAR/rt_getters.h
@@ -104,10 +104,9 @@ __attribute__((always_inline)) INLINE static void rt_get_pressure_tensor(
return;
}
- const float f = rt_params.reduced_speed_of_light_inverse * Fnorm / U[0];
- /* f^2 mustn't be > 1. This may happen since we use the reduced speed of
- * light. */
- const float f2 = min(1.f, f * f);
+ /* f mustn't be > 1. This may happen since we use the reduced speed of light. */
+ const float f = min(1.f, rt_params.reduced_speed_of_light_inverse * Fnorm / U[0]);
+ const float f2 = f * f;
const float rootterm = 4.f - 3.f * f2;
const float chi = (3.f + 4.f * f2) / (5.f + 2.f * sqrtf(rootterm));
diff --git a/src/rt/GEAR/rt_interaction_rates.h b/src/rt/GEAR/rt_interaction_rates.h
index 62db9467a1d32546a7638020f07859db2d341496..fb89bbf6426d31a780bb59fea641264aaae7e703 100644
--- a/src/rt/GEAR/rt_interaction_rates.h
+++ b/src/rt/GEAR/rt_interaction_rates.h
@@ -300,6 +300,7 @@ static void rt_interaction_rates_init(struct rt_props *restrict rt_props,
integration_params.kB,
integration_params.h_planck);
} else {
+ nu_stop_final = -1.;
error("Unknown stellar spectrum type %d", rt_props->stellar_spectrum_type);
}
diff --git a/src/rt/GEAR/rt_slope_limiters_face.h b/src/rt/GEAR/rt_slope_limiters_face.h
index cfde3ea8cbcdc7ddca53cb674aa7541ec4258a65..2b4f8ad939064010554d71f77ae85c12fb528ea0 100644
--- a/src/rt/GEAR/rt_slope_limiters_face.h
+++ b/src/rt/GEAR/rt_slope_limiters_face.h
@@ -126,7 +126,7 @@ __attribute__((always_inline)) INLINE static float rt_limiter_mc(
__attribute__((always_inline)) INLINE static float rt_limiter_vanLeer(
const float dQi, const float dQj) {
const float r = dQj == 0.f ? dQi * 1e6 : dQi / dQj;
- const float absr = fabs(r);
+ const float absr = fabsf(r);
return (r + absr) / (1.f + absr);
}
diff --git a/src/rt/GEAR/rt_thermochemistry.h b/src/rt/GEAR/rt_thermochemistry.h
index 427a64ad94287f8dcfeb554cfe05bba13f492252..d50e9cb62262a599df4c9e7ad1b074257efaa730 100644
--- a/src/rt/GEAR/rt_thermochemistry.h
+++ b/src/rt/GEAR/rt_thermochemistry.h
@@ -189,13 +189,14 @@ static void rt_do_thermochemistry(struct part* restrict p,
/* update radiation fields */
for (int g = 0; g < RT_NGROUPS; g++) {
+ const float e_old = p->rt_data.conserved[g].energy;
const float factor_new = (1.f - dt * rates_by_frequency_bin[g]);
p->rt_data.conserved[g].energy *= factor_new;
for (int i = 0; i < 3; i++) {
p->rt_data.conserved[g].flux[i] *= factor_new;
}
rt_check_unphysical_conserved(&p->rt_data.conserved[g].energy,
- p->rt_data.conserved[g].flux, 2);
+ p->rt_data.conserved[g].flux, e_old, 2);
}
/* copy updated grackle data to particle */
diff --git a/src/rt/GEAR/rt_unphysical.h b/src/rt/GEAR/rt_unphysical.h
index dd1ae1ed33f451db12217e1a31834f91e73a14b0..7804ef16a1216d4e11e0abe446e4085a771e0a33 100644
--- a/src/rt/GEAR/rt_unphysical.h
+++ b/src/rt/GEAR/rt_unphysical.h
@@ -79,16 +79,24 @@ __attribute__((always_inline)) INLINE static void rt_check_unphysical_density(
*
* @param energy pointer to the photon energy
* @param flux pointer to photon fluxes (3 dimensional)
+ * @param e_old photon energy before the change to be checked, if available.
* @param c integer indentifier where this function was called from
*/
__attribute__((always_inline)) INLINE static void rt_check_unphysical_conserved(
- float* energy, float* flux, int c) {
+ float* energy, float* flux, const float e_old, int c) {
/* Check for negative energies */
+ /* Note to self for printouts: Maximal allowable F = E * c.
+ * In some cases, e.g. while cooling, we don't modify the fluxes,
+ * so you can get an estimate of what the photon energy used to be
+ * by dividing the printed out fluxes by the speed of light in
+ * code units */
#ifdef SWIFT_DEBUG_CHECKS
- if (*energy < 0.f && fabs(*energy) > 1.e-1)
- message("Fixing unphysical energy case %d | %.6e | %.6e %.6e %.6e", c,
- *energy, flux[0], flux[1], flux[2]);
+ float ratio = 1.;
+ if (e_old != 0.f) ratio = fabs(*energy / e_old);
+ if (*energy < 0.f && ratio > 1.e-4)
+ message("Fixing unphysical energy case %d | %.6e | %.6e %.6e %.6e | %.6e", c,
+ *energy, flux[0], flux[1], flux[2], ratio);
#endif
if (isinf(*energy) || isnan(*energy))
error("Got inf/nan radiation energy case %d | %.6e | %.6e %.6e %.6e", c,