diff --git a/examples/RadiativeTransferTests/CosmoUniformBox_3D/README b/examples/RadiativeTransferTests/CosmoUniformBox_3D/README
index a5dbcb88d1942dc08660d255d47e8c4813a1292e..0c2a147a53c08de58ef854f148953154920c8e4e 100644
--- a/examples/RadiativeTransferTests/CosmoUniformBox_3D/README
+++ b/examples/RadiativeTransferTests/CosmoUniformBox_3D/README
@@ -2,5 +2,5 @@ A simple box with static particles on a grid, and uniform radiation energy densi
Test that total radiation energy and radiation energy density dilute with the correct scale-factor relation.
The ICs are created to be compatible with GEAR_RT. Recommended configuration:
- --with-rt=GEAR_1 --with-rt-riemann-solver=GLF --with-hydro-dimension=3 --with-hydro=gizmo-mfv \
+ --with-rt=GEAR_11 --with-rt-riemann-solver=GLF --with-hydro-dimension=3 --with-hydro=gizmo-mfv \
--with-riemann-solver=hllc --with-stars=GEAR --with-feedback=none --with-grackle=$GRACKLE_ROOT
diff --git a/examples/RadiativeTransferTests/CosmoUniformBox_3D/makeIC.py b/examples/RadiativeTransferTests/CosmoUniformBox_3D/makeIC.py
index 647080dce6415abcfa37980eb434c296d241c37f..aee6dc2f158488cae735690fddcc7bf2c5db07e3 100755
--- a/examples/RadiativeTransferTests/CosmoUniformBox_3D/makeIC.py
+++ b/examples/RadiativeTransferTests/CosmoUniformBox_3D/makeIC.py
@@ -9,14 +9,14 @@ from swiftsimio.units import cosmo_units
# Unit system we're working with
unitsystem = cosmo_units
-# Box is 100 Mpc in each direction
+# Box is 260 Mpc in each direction
boxsize = 260 * unyt.Mpc
boxsize = boxsize.to(unitsystem["length"])
reduced_speed_of_light_fraction = 1.0
# Number of photon groups
-nPhotonGroups = 1
+nPhotonGroups = 11
# Number of particles in each dimension
# Total number of particles is thus n_p^3
@@ -104,15 +104,16 @@ if __name__ in ("__main__"):
nparts = header.attrs["NumPart_ThisFile"][0]
parts = F["/PartType0"]
- # Generate initial conditions
- E, fluxes = initial_condition(unitsystem)
+ for i in range(nPhotonGroups):
+ # Generate initial conditions
+ E, fluxes = initial_condition(unitsystem)
- # Create photon energy data entry
- dsetname = "PhotonEnergiesGroup1"
- energydata = np.zeros((nparts), dtype=np.float32)
- parts.create_dataset(dsetname, data=E)
+ # Create photon energy data entry
+ dsetname = f"PhotonEnergiesGroup{i+1}"
+ energydata = np.zeros((nparts), dtype=np.float32)
+ parts.create_dataset(dsetname, data=E)
- # Create photon fluxes data entry
- dsetname = "PhotonFluxesGroup1"
- fluxdata = np.zeros((nparts, 3), dtype=np.float32)
- parts.create_dataset(dsetname, data=fluxes)
+ # Create photon fluxes data entry
+ dsetname = f"PhotonFluxesGroup{i+1}"
+ fluxdata = np.zeros((nparts, 3), dtype=np.float32)
+ parts.create_dataset(dsetname, data=fluxes)
diff --git a/examples/RadiativeTransferTests/CosmoUniformBox_3D/plotSolution.py b/examples/RadiativeTransferTests/CosmoUniformBox_3D/plotSolution.py
index d55f837662bf5c5767f5ad1b1e2cbc8298fb61a7..549c16d84313b822be98e2ec9e880163e5988fe6 100755
--- a/examples/RadiativeTransferTests/CosmoUniformBox_3D/plotSolution.py
+++ b/examples/RadiativeTransferTests/CosmoUniformBox_3D/plotSolution.py
@@ -14,6 +14,8 @@ from scipy.optimize import curve_fit as cf
snapshot_base = "output" # snapshot basename
plot_physical_quantities = True
+nPhotonGroups = 11
+
mpl.rcParams["text.usetex"] = True
@@ -65,11 +67,12 @@ def plot_param_over_time(
meta = data.metadata
# Read comoving quantities
- energy = getattr(data.gas.photon_energies, "group1")
+ energy = data.gas.photon_energies.group1
+ for i in range(nPhotonGroups - 1):
+ energy += getattr(data.gas.photon_energies, f"group{i+2}")
mass = data.gas.masses
rho = data.gas.densities
vol = mass / rho
-
energy_density = energy / vol
if plot_physical_quantities:
@@ -83,29 +86,17 @@ def plot_param_over_time(
* np.sum(physical_energy_density)
/ physical_energy_density.shape[0]
)
- analytic_exponent[1] = -3.0
- elif param == "volume":
- plot_param[1].append(1 * np.sum(physical_vol) / physical_vol.shape[0])
- analytic_exponent[1] = 3.0
+ analytic_exponent[1] = -4.0
elif param == "total energy":
plot_param[1].append(1 * np.sum(physical_energy))
- analytic_exponent[1] = 0.0
- elif param == "mass":
- plot_param[1].append(1 * np.sum(physical_mass))
- analytic_exponent[1] = 0.0
+ analytic_exponent[1] = -1.0
if param == "energy density":
plot_param[0].append(1 * np.sum(energy_density) / energy_density.shape[0])
- analytic_exponent[0] = 0.0
- elif param == "volume":
- plot_param[0].append(1 * np.sum(vol) / vol.shape[0])
- analytic_exponent[0] = 0.0
+ analytic_exponent[0] = -1.0
elif param == "total energy":
plot_param[0].append(1 * np.sum(energy))
- analytic_exponent[0] = 0.0
- elif param == "mass":
- plot_param[0].append(1 * np.sum(mass))
- analytic_exponent[0] = 0.0
+ analytic_exponent[0] = -1.0
scale_factor.append(meta.scale_factor)
@@ -117,18 +108,10 @@ def plot_param_over_time(
titles = ["Comoving energy density", "Physical energy density"]
ylabel = "Average energy density"
figname = f"output_energy_density_over_time-{redshift_domain}.png"
- elif param == "volume":
- titles = ["Comoving particle volume", "Physical particle volume"]
- ylabel = "Average particle volume"
- figname = f"output_volume_over_time-{redshift_domain}.png"
elif param == "total energy":
titles = ["Comoving total energy", "Physical total energy"]
ylabel = "Total energy"
figname = f"output_total_energy_over_time-{redshift_domain}.png"
- elif param == "mass":
- titles = ["Comoving total mass", "Physical total mass"]
- ylabel = "Total mass"
- figname = f"output_total_mass_over_time-{redshift_domain}.png"
for i in range(1 + plot_physical_quantities):
ax = fig.add_subplot(1, (1 + plot_physical_quantities), (1 + i))
@@ -147,7 +130,7 @@ def plot_param_over_time(
)
ax.legend()
- ax.set_title(titles[i])
+ ax.set_title(titles[i] + " sum of all groups")
ax.set_xlabel("Scale factor")
secax = ax.secondary_xaxis("top", functions=(a2z, z2a))
@@ -187,10 +170,9 @@ if __name__ in ("__main__"):
sys.exit(1)
snaplist = get_snapshot_list(snapshot_base + f"_{redshift_domain}")
-
if len(snaplist) < 1:
print("No snapshots found!")
exit(1)
- for param in ["energy density", "volume", "total energy", "mass"]:
+ for param in ["energy density", "total energy"]:
plot_param_over_time(snaplist, param, redshift_domain)
diff --git a/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_high_redshift.yml b/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_high_redshift.yml
index f4177d6804a5865a530e5085f8f7032144ea4b98..a258a9d92d5dbc2e70a6b1897e74bbc848e75e59 100644
--- a/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_high_redshift.yml
+++ b/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_high_redshift.yml
@@ -14,8 +14,6 @@ TimeIntegration:
max_nr_rt_subcycles: 1
dt_min: 1.e-17 # The minimal time-step size of the simulation (in internal units).
dt_max: 1.e-2 # The maximal time-step size of the simulation (in internal units).
- time_begin: 0
- time_end: 1.
# Parameters governing the snapshots
Snapshots:
@@ -46,11 +44,12 @@ GEARRT:
f_reduce_c: 1. # reduce the speed of light for the RT solver by multiplying c with this factor
f_limit_cooling_time: 0.0 # (Optional) multiply the cooling time by this factor when estimating maximal next time step. Set to 0.0 to turn computation of cooling time off.
CFL_condition: 0.99 # CFL condition for RT, independent of hydro
- photon_groups_Hz: [1.] # Lower photon frequency group bin edges in Hz. Needs to have exactly N elements, where N is the configured number of bins --with-RT=GEAR_N
+ photon_groups_Hz: [0, 3.2880e+15, 6.5760e+15, 9.8640e+15, 1.3152e+16, 1.6440e+16, 1.9728e+16, 2.3016e+16, 2.6304e+16, 2.9592e+16, 3.2880e+16] # Hz
stellar_luminosity_model: const # Which model to use to determine the stellar luminosities.
- const_stellar_luminosities_LSol: [0.] # (Conditional) constant star luminosities for each photon frequency group to use if stellar_luminosity_model:const is set, in units of Solar Luminosity.
+ const_stellar_luminosities_LSol: [0, 2.221e+04, 2.020e+04, 1.153e+04, 5.122e+03, 1.952e+03, 6.705e+02, 2.140e+02, 6.461e+01, 1.869e+01, 7.158e+00]
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_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.e6 # K
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
diff --git a/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_low_redshift.yml b/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_low_redshift.yml
index 6bac120ae6d9d4e8dcf0239a1e656073d52b97fa..fc07cfe7dd5da62b99f6e54afc116eb99b3cf8a2 100644
--- a/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_low_redshift.yml
+++ b/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_low_redshift.yml
@@ -14,8 +14,6 @@ TimeIntegration:
max_nr_rt_subcycles: 1
dt_min: 1.e-17 # The minimal time-step size of the simulation (in internal units).
dt_max: 1.e-2 # The maximal time-step size of the simulation (in internal units).
- time_begin: 0
- time_end: 1.
# Parameters governing the snapshots
Snapshots:
@@ -46,12 +44,12 @@ GEARRT:
f_reduce_c: 1. # reduce the speed of light for the RT solver by multiplying c with this factor
f_limit_cooling_time: 0.0 # (Optional) multiply the cooling time by this factor when estimating maximal next time step. Set to 0.0 to turn computation of cooling time off.
CFL_condition: 0.99 # CFL condition for RT, independent of hydro
- photon_groups_Hz: [1.] # Lower photon frequency group bin edges in Hz. Needs to have exactly N elements, where N is the configured number of bins --with-RT=GEAR_N
+ photon_groups_Hz: [0, 3.2880e+15, 6.5760e+15, 9.8640e+15, 1.3152e+16, 1.6440e+16, 1.9728e+16, 2.3016e+16, 2.6304e+16, 2.9592e+16, 3.2880e+16] # Hz
stellar_luminosity_model: const # Which model to use to determine the stellar luminosities.
- const_stellar_luminosities_LSol: [0.] # (Conditional) constant star luminosities for each photon frequency group to use if stellar_luminosity_model:const is set, in units of Solar Luminosity.
- hydrogen_mass_fraction: 0.76 # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
- 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.
+ const_stellar_luminosities_LSol: [0, 2.221e+04, 2.020e+04, 1.153e+04, 5.122e+03, 1.952e+03, 6.705e+02, 2.140e+02, 6.461e+01, 1.869e+01, 7.158e+00]
+ hydrogen_mass_fraction: 0.76 # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
+ 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.e6 # (Optinal) Blackbody temperature of the 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
diff --git a/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_medium_redshift.yml b/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_medium_redshift.yml
index 9189a2ca73dbfa70d48a3e9b1765b82c8ab8e864..412be5b56e32109ec11def92dc349f1591d57581 100644
--- a/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_medium_redshift.yml
+++ b/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_medium_redshift.yml
@@ -14,8 +14,6 @@ TimeIntegration:
max_nr_rt_subcycles: 1
dt_min: 1.e-17 # The minimal time-step size of the simulation (in internal units).
dt_max: 1.e-2 # The maximal time-step size of the simulation (in internal units).
- time_begin: 0
- time_end: 1.
# Parameters governing the snapshots
Snapshots:
@@ -46,12 +44,12 @@ GEARRT:
f_reduce_c: 1. # reduce the speed of light for the RT solver by multiplying c with this factor
f_limit_cooling_time: 0.0 # (Optional) multiply the cooling time by this factor when estimating maximal next time step. Set to 0.0 to turn computation of cooling time off.
CFL_condition: 0.99 # CFL condition for RT, independent of hydro
- photon_groups_Hz: [1.] # Lower photon frequency group bin edges in Hz. Needs to have exactly N elements, where N is the configured number of bins --with-RT=GEAR_N
+ photon_groups_Hz: [0, 3.2880e+15, 6.5760e+15, 9.8640e+15, 1.3152e+16, 1.6440e+16, 1.9728e+16, 2.3016e+16, 2.6304e+16, 2.9592e+16, 3.2880e+16] # Hz
stellar_luminosity_model: const # Which model to use to determine the stellar luminosities.
- const_stellar_luminosities_LSol: [0.] # (Conditional) constant star luminosities for each photon frequency group to use if stellar_luminosity_model:const is set, in units of Solar Luminosity.
- hydrogen_mass_fraction: 0.76 # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
- 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.
+ const_stellar_luminosities_LSol: [0, 2.221e+04, 2.020e+04, 1.153e+04, 5.122e+03, 1.952e+03, 6.705e+02, 2.140e+02, 6.461e+01, 1.869e+01, 7.158e+00]
+ hydrogen_mass_fraction: 0.76 # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
+ 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.e6 # (Optinal) Blackbody temperature of the 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
diff --git a/src/rt/GEAR/rt.h b/src/rt/GEAR/rt.h
index f2c2f31e7f6f709c591fad8fdb62494b7776e6dc..802a1fa449fec69b54580c84d5ad87dd8658dea1 100644
--- a/src/rt/GEAR/rt.h
+++ b/src/rt/GEAR/rt.h
@@ -462,7 +462,7 @@ __attribute__((always_inline)) INLINE static void rt_end_gradient(
* @param cosmo #cosmology data structure.
*/
__attribute__((always_inline)) INLINE static void rt_finalise_transport(
- struct part* restrict p, const double dt,
+ struct part* restrict p, struct rt_props* rtp, const double dt,
const struct cosmology* restrict cosmo) {
#ifdef SWIFT_RT_DEBUG_CHECKS
@@ -480,18 +480,35 @@ __attribute__((always_inline)) INLINE static void rt_finalise_transport(
struct rt_part_data* restrict rtd = &p->rt_data;
const float Vinv = 1.f / p->geometry.volume;
+ /* Do not redshift if we have a constant spectrum (type == 0) */
+ const float redshift_factor =
+ (rtp->stellar_spectrum_type == 0) ? 0. : cosmo->H * dt;
+
for (int g = 0; g < RT_NGROUPS; g++) {
const float e_old = rtd->radiation[g].energy_density;
/* Note: in this scheme, we're updating d/dt (U * V) + sum F * A * dt = 0.
* So we'll need the division by the volume here. */
+
rtd->radiation[g].energy_density += rtd->flux[g].energy * Vinv;
+ rtd->radiation[g].energy_density -=
+ rtd->radiation[g].energy_density *
+ redshift_factor; // Energy lost due to redshift
rtd->radiation[g].flux[0] += rtd->flux[g].flux[0] * Vinv;
+ rtd->radiation[g].flux[0] -=
+ rtd->radiation[g].flux[0] *
+ redshift_factor; // Energy lost due to redshift
rtd->radiation[g].flux[1] += rtd->flux[g].flux[1] * Vinv;
+ rtd->radiation[g].flux[1] -=
+ rtd->radiation[g].flux[1] *
+ redshift_factor; // Energy lost due to redshift
rtd->radiation[g].flux[2] += rtd->flux[g].flux[2] * Vinv;
+ rtd->radiation[g].flux[2] -=
+ rtd->radiation[g].flux[2] *
+ redshift_factor; // Energy lost due to redshift
rt_check_unphysical_state(&rtd->radiation[g].energy_density,
rtd->radiation[g].flux, e_old, /*callloc=*/4);
diff --git a/src/rt/SPHM1RT/rt.h b/src/rt/SPHM1RT/rt.h
index 64ad8eaa263d20c609352ec48ca293f2e3be3bf1..0cb2b94202158d5f63df771c71e2217e08dd439d 100644
--- a/src/rt/SPHM1RT/rt.h
+++ b/src/rt/SPHM1RT/rt.h
@@ -451,7 +451,7 @@ __attribute__((always_inline)) INLINE static void rt_end_gradient(
* @param cosmo #cosmology data structure.
*/
__attribute__((always_inline)) INLINE static void rt_finalise_transport(
- struct part* restrict p, const double dt,
+ struct part* restrict p, struct rt_props* rtp, const double dt,
const struct cosmology* restrict cosmo) {
struct rt_part_data* rpd = &p->rt_data;
diff --git a/src/rt/debug/rt.h b/src/rt/debug/rt.h
index 3742ec1129dc348dfe27160e28c06e3d84236ebf..cf14510431988fa127d1fbff06c0d066ea9099b9 100644
--- a/src/rt/debug/rt.h
+++ b/src/rt/debug/rt.h
@@ -300,9 +300,8 @@ __attribute__((always_inline)) INLINE static void rt_end_gradient(
* @param cosmo #cosmology data structure.
*/
__attribute__((always_inline)) INLINE static void rt_finalise_transport(
- struct part* restrict p, const double dt,
+ struct part* restrict p, struct rt_props* rtp, const double dt,
const struct cosmology* restrict cosmo) {
-
rt_debug_sequence_check(p, 3, __func__);
if (p->rt_data.debug_calls_iact_transport_interaction == 0)
diff --git a/src/rt/none/rt.h b/src/rt/none/rt.h
index 9edc956e1885080ac98f526b1f6f315c8267694a..6254cb4040b8d530141e5bd6ad0e9d169f3229ad 100644
--- a/src/rt/none/rt.h
+++ b/src/rt/none/rt.h
@@ -248,9 +248,8 @@ __attribute__((always_inline)) INLINE static void rt_end_gradient(
* @param cosmo #cosmology data structure.
*/
__attribute__((always_inline)) INLINE static void rt_finalise_transport(
- struct part* restrict p, const double dt,
+ struct part* restrict p, struct rt_props* rtp, const double dt,
const struct cosmology* restrict cosmo) {}
-
/**
* @brief Do the thermochemistry on a particle.
*
diff --git a/src/runner_others.c b/src/runner_others.c
index 65da2188da78f77440493b9715f74cf788f96735..d7b79c9870415409e51911e719f890fc0d3bfe67 100644
--- a/src/runner_others.c
+++ b/src/runner_others.c
@@ -1144,7 +1144,7 @@ void runner_do_rt_tchem(struct runner *r, struct cell *c, int timer) {
error("Got part with negative time-step: %lld, %.6g", p->id, dt);
#endif
- rt_finalise_transport(p, dt, cosmo);
+ rt_finalise_transport(p, rt_props, dt, cosmo);
/* And finally do thermochemistry */
rt_tchem(p, xp, rt_props, cosmo, hydro_props, phys_const, us, dt);