diff --git a/configure.ac b/configure.ac
index a7df903a32f8aa04a6eee9c823e56e4a9df52fee..faa1626a69a23f0d6803af13179793ff3540e116 100644
--- a/configure.ac
+++ b/configure.ac
@@ -2583,26 +2583,6 @@ case "$with_rt" in
esac
-AC_ARG_ENABLE([rt-hydro-controlled-injection],
- [AS_HELP_STRING([--enable-rt-hydro-controlled-injection],
- [enable the radiative transfer injection scheme where \
- active hydro particles gather radiation from neighbouring \
- stars. Experimental feature, use with caution. \
- @<:@no/yes@:>@ Default: no]
- )],
- [with_hydro_controlled_injection="${enableval}"],
- [with_hydro_controlled_injection="no"]
-)
-
-rt_extra_msg="" # extra message for ./configure printout
-if test "$with_hydro_controlled_injection" = "yes"; then
- rt_extra_msg=", hydro controlled injection"
- AC_DEFINE([RT_HYDRO_CONTROLLED_INJECTION], 1,
- [Enable hydro controlled radiative transfer injection scheme])
-fi
-
-
-
# Check for git, needed for revision stamps.
AC_PATH_PROG([GIT_CMD], [git])
AC_SUBST([GIT_CMD])
@@ -2745,7 +2725,7 @@ AC_MSG_RESULT([
Star feedback model : $with_feedback_name
Sink particle model : $with_sink
Black holes model : $with_black_holes
- Radiative transfer : $with_rt$rt_extra_msg
+ Radiative transfer : $with_rt
Atomic operations in tasks : $enable_atomics_within_tasks
Individual timers : $enable_timers
diff --git a/examples/RadiativeTransferTests/Advection_2D/rt_advection2D.yml b/examples/RadiativeTransferTests/Advection_2D/rt_advection2D.yml
index 28b1483c851e3a3cf2e14ecb2fe5b4e5f16c4e85..c143bf77d9e420cfe690ebd06c7ec57a3c1bfe12 100644
--- a/examples/RadiativeTransferTests/Advection_2D/rt_advection2D.yml
+++ b/examples/RadiativeTransferTests/Advection_2D/rt_advection2D.yml
@@ -39,7 +39,7 @@ InitialConditions:
periodic: 1 # peridoc ICs
Scheduler:
- max_top_level_cells: 128
+ max_top_level_cells: 32
# Parameters for the self-gravity scheme
Gravity:
diff --git a/examples/RadiativeTransferTests/RandomizedBox_3D/randomized-rt.yml b/examples/RadiativeTransferTests/RandomizedBox_3D/randomized-rt.yml
index 7c67e8b5e9635bd6faf0980cc4a107a862394086..47194c52626db1de181ad7c919d9cee0a790a3bc 100644
--- a/examples/RadiativeTransferTests/RandomizedBox_3D/randomized-rt.yml
+++ b/examples/RadiativeTransferTests/RandomizedBox_3D/randomized-rt.yml
@@ -45,16 +45,19 @@ Stars:
Scheduler:
max_top_level_cells: 64
+ dependency_graph_frequency: 0 # (Optional) Dumping frequency of the dependency graph. By default, writes only at the first step.
Restarts:
delta_hours: 72 # (Optional) decimal hours between dumps of restart files.
+ enable: 1 # (Optional) whether to enable dumping restarts at fixed intervals.
+ stop_steps: 500 # (Optional) how many steps to process before checking if the <subdir>/stop file exists. When present the application will attempt to exit early, dumping restart files first.
# Parameters for the self-gravity scheme
Gravity:
- mesh_side_length: 12 # Number of cells along each axis for the periodic gravity mesh.
+ mesh_side_length: 32 # Number of cells along each axis for the periodic gravity mesh.
eta: 0.025 # Constant dimensionless multiplier for time integration.
MAC: adaptive # Choice of mulitpole acceptance criterion: 'adaptive' OR 'geometric'.
- epsilon_fmm: 0.01 # Tolerance parameter for the adaptive multipole acceptance criterion.
+ epsilon_fmm: 0.001 # Tolerance parameter for the adaptive multipole acceptance criterion.
theta_cr: 0.7 # Opening angle for the purely gemoetric criterion.
# comoving_DM_softening: 0.0026994 # Comoving Plummer-equivalent softening length for DM particles (in internal units).
# max_physical_DM_softening: 0.0007 # Maximal Plummer-equivalent softening length in physical coordinates for DM particles (in internal units).
diff --git a/examples/RadiativeTransferTests/RandomizedBox_3D/run.sh b/examples/RadiativeTransferTests/RandomizedBox_3D/run.sh
index 7ab69b15f5a35c0aee1fca4cd3ede2bf07bc3745..5446a4c5b212461c1d4a4626ee624d03f5e430e4 100755
--- a/examples/RadiativeTransferTests/RandomizedBox_3D/run.sh
+++ b/examples/RadiativeTransferTests/RandomizedBox_3D/run.sh
@@ -9,8 +9,20 @@ if [ ! -f 'randomized-sine.hdf5' ]; then
python3 makeIC.py
fi
+cmd=../../swift
+if [ $# -gt 0 ]; then
+ case "$1" in
+ g | gdb)
+ cmd='gdb --args ../../swift'
+ ;;
+ *)
+ echo unknown cmdline param, running without gdb
+ ;;
+ esac
+fi
+
# Run SWIFT with RT
-../../swift \
+$cmd \
--hydro \
--threads=9 \
--verbose=0 \
diff --git a/examples/RadiativeTransferTests/UniformBox_3D/rt_sanity_checks-GEAR.py b/examples/RadiativeTransferTests/UniformBox_3D/rt_sanity_checks-GEAR.py
index 051b5db1d743e1b112fca3072752dfdee1b315ad..1f003bd0e54e255962b292b1fe4504fa698f06ec 100755
--- a/examples/RadiativeTransferTests/UniformBox_3D/rt_sanity_checks-GEAR.py
+++ b/examples/RadiativeTransferTests/UniformBox_3D/rt_sanity_checks-GEAR.py
@@ -244,7 +244,7 @@ def check_injection(snapdata, rundata):
np.zeros(rundata.ngroups, dtype=np.float64) * unyt.erg
)
- if rundata.use_const_emission_rate and not rundata.hydro_controlled_injection:
+ if rundata.use_const_emission_rate:
if len(snapdata) <= 2:
# because it's useless to check only snap_0000
print("Check 1b: You need at least 2 snapshots to do this particular test")
@@ -324,7 +324,7 @@ def check_injection(snapdata, rundata):
# Create additional plots?
# --------------------------------
- if rundata.use_const_emission_rate and not rundata.hydro_controlled_injection:
+ if rundata.use_const_emission_rate:
if not skip_plots and len(energies_for_plot) > 2:
# Show me the plot that the injected energy
# is correctly bounded
diff --git a/examples/RadiativeTransferTests/UniformBox_3D/rt_sanity_checks.py b/examples/RadiativeTransferTests/UniformBox_3D/rt_sanity_checks.py
index 86eafab768de3b1482bc6fdbb6fd4b357a6f6981..f28dcc661d5d7dad9b97cfc6c89f60658c1fae1d 100755
--- a/examples/RadiativeTransferTests/UniformBox_3D/rt_sanity_checks.py
+++ b/examples/RadiativeTransferTests/UniformBox_3D/rt_sanity_checks.py
@@ -49,8 +49,6 @@ skip_last_snap = False # skip snap_0max.hdf5
print_diffs = True # print differences you find
break_on_diff = False # quit when you find a difference
-hydro_controlled_injection = False
-
if len(argv) > 1:
file_prefix = argv[1]
@@ -230,21 +228,24 @@ def check_stars_sanity(snapdata):
this = snap.stars
nspart = snapdata[0].stars.coords.shape[0]
- if hydro_controlled_injection:
- if (this.EmissionRateSet != 1).any():
- print("- checking stars sanity pt2", snap.snapnr)
- print("--- Emisison Rates not consistent")
- count = 0
- for i in range(nspart):
- if this.EmissionRateSet[i] != 1:
- count += 1
- if print_diffs:
- print("-----", this.EmissionRateSet[i], "ID", this.IDs[i])
-
- print("--- count", count, "/", this.EmissionRateSet.shape[0])
+ fishy = np.logical_and(this.EmissionRateSet == 0, this.RadiationEmittedTot > 0)
+ if fishy.any():
+ print("- checking stars sanity pt1", snap.snapnr)
+ print("--- Emitted with unset emission rate")
+ for i in range(nspart):
+ if this.EmissionRateSet[i] != 1:
+ count += 1
+ if print_diffs:
+ print("----- IDs", this.IDs[fishy])
+ print("----- EmissionRateSet", this.EmissionRateSet[fishy])
+ print(
+ "----- RadiationEmittedTot", this.RadiationEmittedTot[fishy]
+ )
+
+ print("--- count", count, "/", this.EmissionRateSet.shape[0])
- if break_on_diff:
- quit()
+ if break_on_diff:
+ quit()
return
@@ -274,6 +275,7 @@ def check_stars_hydro_interaction_sanity(snapdata):
if snap.has_stars:
sum_star_tot_radiation = stars.RadiationEmittedTot.sum()
else:
+ print("Found no stars")
sum_star_tot_radiation = 0.0
if sum_gas_tot_radiation != sum_star_tot_radiation:
@@ -292,25 +294,28 @@ def check_stars_hydro_interaction_sanity(snapdata):
# --------------------------------------------------------------
# check that we have the correct amount of interactions recorded
# for injection prep between stars and gas
- # --------------------------------------------------------------
- sum_gas_tot_prep = gas.InjectPrepCountsTot.sum()
- if snap.has_stars:
- sum_star_tot_prep = stars.InjectPrepCountsTot.sum()
- else:
- sum_star_tot_prep = 0.0
-
- if sum_gas_tot_prep != sum_star_tot_prep:
- print("- checking hydro v star sanity pt2; snapshot", snap.snapnr)
- print(
- "--- Total interactions between gas and stars in prep is wrong:",
- sum_gas_tot_prep,
- "stars",
- sum_star_tot_prep,
- "diff",
- sum_star_tot_prep - sum_gas_tot_prep,
- )
- if break_on_diff:
- quit()
+ # !! Can't do this check any longer since we moved the injection
+ # !! prep to the star density loop. The repeats and resets in
+ # !! the ghost mess the counts between parts and sparts up.
+ # ---------------------------------------------------------------
+ # sum_gas_tot_prep = gas.InjectPrepCountsTot.sum()
+ # if snap.has_stars:
+ # sum_star_tot_prep = stars.InjectPrepCountsTot.sum()
+ # else:
+ # sum_star_tot_prep = 0.0
+ #
+ # if sum_gas_tot_prep != sum_star_tot_prep:
+ # print("- checking hydro v star sanity pt2; snapshot", snap.snapnr)
+ # print(
+ # "--- Total interactions between gas and stars in prep is wrong:",
+ # sum_gas_tot_prep,
+ # "stars",
+ # sum_star_tot_prep,
+ # "diff",
+ # sum_star_tot_prep - sum_gas_tot_prep,
+ # )
+ # if break_on_diff:
+ # quit()
return
@@ -323,8 +328,6 @@ def main():
snapdata, rundata = get_snap_data(
prefix=file_prefix, skip_snap_zero=skip_snap_zero, skip_last_snap=skip_last_snap
)
- global hydro_controlled_injection
- hydro_controlled_injection = rundata.hydro_controlled_injection
check_hydro_sanity(snapdata)
check_stars_sanity(snapdata)
diff --git a/examples/RadiativeTransferTests/UniformBox_3D/rt_uniform_box_checks.py b/examples/RadiativeTransferTests/UniformBox_3D/rt_uniform_box_checks.py
index e47fae3ef394942fcadfb882db72dbb5fc039429..c8a16ed9ed82d05d337ff000881a98149efa155f 100755
--- a/examples/RadiativeTransferTests/UniformBox_3D/rt_uniform_box_checks.py
+++ b/examples/RadiativeTransferTests/UniformBox_3D/rt_uniform_box_checks.py
@@ -262,18 +262,59 @@ def check_all_stars_is_equal(snapdata):
# Smoothing Lengths
if not skip_sml:
- diff = np.abs((ref.gas.h - compare.gas.h) / ref.gas.h)
+ diff = np.abs((ref.stars.h - compare.stars.h) / ref.stars.h)
if (diff > float_comparison_tolerance).any():
print("- Comparing stars", ref.snapnr, "->", compare.snapnr)
print("--- Smoothing Lengths vary")
if print_diffs:
for i in range(npart):
- if ((ref.gas.h[i] - compare.gas.h[i]) / ref.gas.h[i]).any():
- print(ref.gas.h[i], "|", compare.gas.h[i])
+ if (
+ (ref.stars.h[i] - compare.stars.h[i]) / ref.stars.h[i]
+ ).any():
+ print(ref.stars.h[i], "|", compare.stars.h[i])
if break_on_diff:
quit()
+ # Check all emission rates are set everywhere
+ fishy = ref.stars.EmissionRateSet != compare.stars.EmissionRateSet
+ if fishy.any():
+
+ print("- Comparing stars", ref.snapnr, "->", compare.snapnr)
+ print("--- EmissionRateSet vary")
+ if print_diffs:
+ for i in range(npart):
+ if ref.stars.EmissionRateSet[i] != compare.stars.EmissionRateSet[i]:
+ print(
+ ref.stars.EmissionRateSet[i],
+ "|",
+ compare.stars.EmissionRateSet[i],
+ )
+
+ if break_on_diff:
+ quit()
+
+ # Check all emitted radiation is equal
+ fishy = ref.stars.InjectionInteractions != compare.stars.InjectionInteractions
+ if fishy.any():
+
+ print("- Comparing stars", ref.snapnr, "->", compare.snapnr)
+ print("--- InjectionInteractions vary")
+ if print_diffs:
+ for i in range(npart):
+ if (
+ ref.stars.InjectionInteractions[i]
+ != compare.stars.InjectionInteractions[i]
+ ):
+ print(
+ ref.stars.InjectionInteractions[i],
+ "|",
+ compare.stars.InjectionInteractions[i],
+ )
+
+ if break_on_diff:
+ quit()
+
return
diff --git a/examples/RadiativeTransferTests/UniformBox_3D/run.sh b/examples/RadiativeTransferTests/UniformBox_3D/run.sh
index 21762615467eb9457e7361ed2ef736177dd0c1b7..e3cebd4932ef1b28bca6077d0cdcc051d2c83875 100755
--- a/examples/RadiativeTransferTests/UniformBox_3D/run.sh
+++ b/examples/RadiativeTransferTests/UniformBox_3D/run.sh
@@ -9,8 +9,20 @@ if [ ! -f 'uniformBox-rt.hdf5' ]; then
python3 makeIC.py
fi
+cmd=../../swift
+if [ $# -gt 0 ]; then
+ case "$1" in
+ g | gdb)
+ cmd='gdb --args ../../swift'
+ ;;
+ *)
+ echo unknown cmdline param, running without gdb
+ ;;
+ esac
+fi
+
# Run SWIFT with RT
-../../swift \
+$cmd \
--hydro --threads=4 --stars --external-gravity \
--feedback --radiation \
uniform_rt_timestep_output_sync.yml 2>&1 | tee output.log
diff --git a/examples/RadiativeTransferTests/UniformBox_3D/swift_rt_GEAR_io.py b/examples/RadiativeTransferTests/UniformBox_3D/swift_rt_GEAR_io.py
index 3e350665bf739cd84ece0ab4f2474383a247cced..a9b427dc19a0bbd7eefa3a8a5ca5b079de9ac64a 100644
--- a/examples/RadiativeTransferTests/UniformBox_3D/swift_rt_GEAR_io.py
+++ b/examples/RadiativeTransferTests/UniformBox_3D/swift_rt_GEAR_io.py
@@ -92,7 +92,6 @@ class Rundata(object):
def __init__(self):
self.units = None
- self.hydro_controlled_injection = False
self.use_const_emission_rate = False
self.has_stars = False # assume we don't have stars, check while reading in
@@ -203,9 +202,6 @@ def get_snap_data(prefix="output", skip_snap_zero=False, skip_last_snap=False):
else:
quit()
- if "hydro controlled" in scheme:
- rundata.hydro_controlled_injection = True
-
rundata.reduced_speed_of_light = firstfile.metadata.reduced_lightspeed
# -------------------
diff --git a/examples/RadiativeTransferTests/UniformBox_3D/swift_rt_debug_io.py b/examples/RadiativeTransferTests/UniformBox_3D/swift_rt_debug_io.py
index 6f3be1a53a083bfc6fe15234719e406cd977a55b..c68f9bd4dce8d9ba29a753ab381bd0b406f7168f 100644
--- a/examples/RadiativeTransferTests/UniformBox_3D/swift_rt_debug_io.py
+++ b/examples/RadiativeTransferTests/UniformBox_3D/swift_rt_debug_io.py
@@ -67,8 +67,8 @@ class RTStarData(object):
self.h = None
self.EmissionRateSet = None
+ self.InjectionInteractions = None
self.RadiationEmittedTot = None
- self.InjectPrepCountsTot = None
return
@@ -94,7 +94,6 @@ class Rundata(object):
"""
def __init__(self):
- self.hydro_controlled_injection = False
self.has_stars = False # assume we don't have stars, check while reading in
return
@@ -167,8 +166,6 @@ def get_snap_data(prefix="output", skip_snap_zero=False, skip_last_snap=False):
"Compile swift --with-rt=debug",
)
- if "hydro controlled" in scheme:
- rundata.hydro_controlled_injection = True
F.close()
for f in hdf5files:
@@ -206,13 +203,17 @@ def get_snap_data(prefix="output", skip_snap_zero=False, skip_last_snap=False):
newsnap.gas.ThermochemistryDone = Gas["RTDebugThermochemistryDone"][:][inds]
newsnap.gas.RadiationAbsorbedTot = Gas["RTDebugRadAbsorbedTot"][:][inds]
- newsnap.gas.InjectPrepCountsTot = Gas["RTDebugStarsInjectPrepTotCounts"][:][
- inds
- ]
+ has_stars = False
try:
Stars = F["PartType4"]
ids = Stars["ParticleIDs"][:]
+ has_stars = True
+ except KeyError:
+ has_stars = False
+
+ if has_stars:
+ newsnap.has_stars = has_stars
inds = np.argsort(ids)
newsnap.stars.IDs = ids[inds]
@@ -220,13 +221,8 @@ def get_snap_data(prefix="output", skip_snap_zero=False, skip_last_snap=False):
newsnap.stars.h = Stars["SmoothingLengths"][:][inds]
newsnap.stars.EmissionRateSet = Stars["RTDebugEmissionRateSet"][:][inds]
-
+ newsnap.stars.InjectionInteractions = Stars["RTDebugHydroIact"][:][inds]
newsnap.stars.RadiationEmittedTot = Stars["RTDebugRadEmittedTot"][:][inds]
- newsnap.stars.InjectPrepCountsTot = Stars[
- "RTDebugHydroInjectPrepCountsTot"
- ][:][inds]
- except KeyError:
- newsnap.has_stars = False
snapdata.append(newsnap)
diff --git a/examples/RadiativeTransferTests/UniformBox_3D/uniform_rt_timestep_output_sync.yml b/examples/RadiativeTransferTests/UniformBox_3D/uniform_rt_timestep_output_sync.yml
index 324de5ca94f7ed5648d6b4669c3604c7c2a17ac9..114b93a8b2f63b4d8c713314510f3cc369a662c7 100644
--- a/examples/RadiativeTransferTests/UniformBox_3D/uniform_rt_timestep_output_sync.yml
+++ b/examples/RadiativeTransferTests/UniformBox_3D/uniform_rt_timestep_output_sync.yml
@@ -25,7 +25,7 @@ Snapshots:
# Parameters governing the conserved quantities statistics
Statistics:
time_first: 0.
- delta_time: 5e-3 # Time between statistics output
+ delta_time: 2.980232e-08 # Time between statistics output
# Parameters for the hydrodynamics scheme
SPH:
@@ -43,10 +43,6 @@ Scheduler:
cell_split_size: 25 # Lower than default to test going deep(er) in the tree
dependency_graph_frequency: 0
-# Debugging/Development RT scheme
-DebugRT:
- all_parts_have_stars: 0 # Set to 1 to do additional tests that only work if all hydro particles have a star particle neighbour
-
GEARRT:
f_reduce_c: 1e-3 # reduce the speed of light for the RT solver by multiplying c with this factor
CFL_condition: 0.99 # CFL condition for RT, independent of hydro
diff --git a/examples/main.c b/examples/main.c
index 7ea08228e1527860ae6f058287cc21a721c57541..a4909ffc39036b39f9e87fd5d86cda6edaef4681 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -1813,7 +1813,7 @@ int main(int argc, char *argv[]) {
if (with_stars) stars_props_clean(e.stars_properties);
if (with_cooling || with_temperature) cooling_clean(e.cooling_func);
if (with_feedback) feedback_clean(e.feedback_props);
- if (with_rt) rt_clean(e.rt_props);
+ if (with_rt) rt_clean(e.rt_props, restart);
engine_clean(&e, /*fof=*/0, restart);
free(params);
free(output_options);
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index c1b474020a880d8ef9e00e017282880ba699b6f3..8f7045d0e0ca46aca63794d2ad2d11159cee9ac0 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -679,10 +679,6 @@ GEARSink:
# Parameters related to radiative transfer ---------------------------------------
-# Debugging/Development RT scheme
-DebugRT:
- all_parts_have_stars: 0 # Set to 1 to do additional tests that only work if all hydro particles have a star particle neighbour
-
GEARRT:
f_reduce_c: 1e-3 # reduce the speed of light for the RT solver by multiplying c with this factor
CFL_condition: 0.9 # CFL condition for RT, independent of hydro
diff --git a/src/Makefile.am b/src/Makefile.am
index f7a385534da87d728e18bce80a531e377fda6a24..eef7fc07d53fe73866139d9009f526331abd7e7a 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -123,7 +123,7 @@ AM_SOURCES += runner_doiact_stars.c runner_doiact_black_holes.c runner_ghost.c
AM_SOURCES += runner_recv.c runner_pack.c
AM_SOURCES += runner_sort.c runner_drift.c runner_black_holes.c runner_time_integration.c
AM_SOURCES += runner_doiact_hydro_vec.c runner_others.c runner_doiact_sinks.c
-AM_SOURCES += runner_doiact_rt.c runner_doiact_sinks_merger.c
+AM_SOURCES += runner_doiact_sinks_merger.c
AM_SOURCES += cell.c cell_convert_part.c cell_drift.c cell_lock.c cell_pack.c cell_split.c
AM_SOURCES += cell_unskip.c
AM_SOURCES += engine.c engine_maketasks.c engine_split_particles.c engine_strays.c
@@ -160,7 +160,7 @@ nobase_noinst_HEADERS += gravity_iact.h kernel_long_gravity.h vector.h accumulat
nobase_noinst_HEADERS += runner_doiact_nosort.h runner_doiact_hydro.h runner_doiact_stars.h runner_doiact_black_holes.h runner_doiact_grav.h
nobase_noinst_HEADERS += runner_doiact_functions_hydro.h runner_doiact_functions_stars.h runner_doiact_functions_black_holes.h
nobase_noinst_HEADERS += runner_doiact_functions_limiter.h runner_doiact_limiter.h units.h intrinsics.h minmax.h
-nobase_noinst_HEADERS += runner_doiact_rt.h runner_doiact_functions_rt.h runner_doiact_sinks.h runner_doiact_functions_sinks.h
+nobase_noinst_HEADERS += runner_doiact_sinks.h runner_doiact_functions_sinks.h
nobase_noinst_HEADERS += runner_doiact_sinks_merger.h runner_doiact_functions_sinks_merger.h
nobase_noinst_HEADERS += kick.h timestep.h drift.h adiabatic_index.h io_properties.h dimension.h part_type.h periodic.h memswap.h
nobase_noinst_HEADERS += timestep_limiter.h timestep_limiter_iact.h timestep_sync.h timestep_sync_part.h timestep_limiter_struct.h
@@ -248,7 +248,6 @@ nobase_noinst_HEADERS += riemann/riemann_hllc.h riemann/riemann_trrs.h
nobase_noinst_HEADERS += riemann/riemann_exact.h riemann/riemann_vacuum.h
nobase_noinst_HEADERS += riemann/riemann_checks.h
nobase_noinst_HEADERS += rt.h
-nobase_noinst_HEADERS += rt_active.h
nobase_noinst_HEADERS += rt_additions.h
nobase_noinst_HEADERS += rt_io.h
nobase_noinst_HEADERS += rt_parameters.h
@@ -269,9 +268,7 @@ nobase_noinst_HEADERS += rt/debug/rt_iact.h
nobase_noinst_HEADERS += rt/debug/rt_io.h
nobase_noinst_HEADERS += rt/debug/rt_parameters.h
nobase_noinst_HEADERS += rt/debug/rt_properties.h
-nobase_noinst_HEADERS += rt/debug/rt_stellar_emission_rate.h
nobase_noinst_HEADERS += rt/debug/rt_struct.h
-nobase_noinst_HEADERS += rt/debug/rt_thermochemistry.h
nobase_noinst_HEADERS += rt/GEAR/rt.h
nobase_noinst_HEADERS += rt/GEAR/rt_additions.h
nobase_noinst_HEADERS += rt/GEAR/rt_blackbody.h
diff --git a/src/cell.c b/src/cell.c
index cf439427d8a70c2b7f29594a3911424902325ca6..9c735106420a33f66ff142f5f3976c75737f6cda 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -579,7 +579,6 @@ void cell_clean_links(struct cell *c, void *data) {
c->hydro.gradient = NULL;
c->hydro.force = NULL;
c->hydro.limiter = NULL;
- c->hydro.rt_inject = NULL;
c->hydro.rt_gradient = NULL;
c->hydro.rt_transport = NULL;
c->grav.grav = NULL;
diff --git a/src/cell.h b/src/cell.h
index 556e7671c56164761e89069e0ba63410b0195e09..fede406a0e6caf5379abce55e6f996ab6ca2cb30 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -574,8 +574,6 @@ void cell_activate_subcell_black_holes_tasks(struct cell *ci, struct cell *cj,
const int with_timestep_sync);
void cell_activate_subcell_external_grav_tasks(struct cell *ci,
struct scheduler *s);
-void cell_activate_subcell_rt_tasks(struct cell *ci, struct cell *cj,
- struct scheduler *s);
void cell_activate_super_spart_drifts(struct cell *c, struct scheduler *s);
void cell_activate_super_sink_drifts(struct cell *c, struct scheduler *s);
void cell_activate_drift_part(struct cell *c, struct scheduler *s);
diff --git a/src/cell_hydro.h b/src/cell_hydro.h
index 26c507d28a56ada1c451cf48082a6b96302f0802..eb24ce84751faeb047aacf0f5ec8c59bf5756c9a 100644
--- a/src/cell_hydro.h
+++ b/src/cell_hydro.h
@@ -112,9 +112,6 @@ struct cell_hydro {
/*! Radiative transfer ghost in task */
struct task *rt_in;
- /*! Task for self/pair injection step of radiative transfer */
- struct link *rt_inject;
-
/*! Radiative transfer ghost1 task (finishes up injection) */
struct task *rt_ghost1;
diff --git a/src/cell_unskip.c b/src/cell_unskip.c
index c4bb8cd916d0f56895365281c7a5067d06b6bc7e..d26f1db2e91500bca4e180917372d088e2c8a371 100644
--- a/src/cell_unskip.c
+++ b/src/cell_unskip.c
@@ -29,7 +29,6 @@
#include "active.h"
#include "engine.h"
#include "feedback.h"
-#include "rt_active.h"
#include "space_getsid.h"
extern int engine_star_resort_task_depth;
@@ -1431,126 +1430,6 @@ void cell_activate_subcell_external_grav_tasks(struct cell *ci,
}
}
-/**
- * @brief Traverse a sub-cell task and activate the stars drift tasks that are
- * required by a Radiative Transfer (injection) task
- *
- * @param ci The first #cell we recurse in.
- * @param cj The second #cell we recurse in.
- * @param s The task #scheduler.
- */
-void cell_activate_subcell_rt_tasks(struct cell *ci, struct cell *cj,
- struct scheduler *s) {
- const struct engine *e = s->space->e;
-
- /* Store the current dx_max and h_max values. */
- ci->stars.dx_max_part_old = ci->stars.dx_max_part;
- ci->stars.h_max_old = ci->stars.h_max;
- ci->hydro.dx_max_part_old = ci->hydro.dx_max_part;
- ci->hydro.h_max_old = ci->hydro.h_max;
-
- if (cj != NULL) {
- cj->stars.dx_max_part_old = cj->stars.dx_max_part;
- cj->stars.h_max_old = cj->stars.h_max;
- cj->hydro.dx_max_part_old = cj->hydro.dx_max_part;
- cj->hydro.h_max_old = cj->hydro.h_max;
- }
-
- /* Self interaction? */
- if (cj == NULL) {
- const int ci_active = rt_should_iact_cell(ci, e);
-
- /* Do anything? */
- if (!ci_active || ci->hydro.count == 0) return;
-
- /* Recurse? */
- if (cell_can_recurse_in_self_stars_task(ci)) {
- /* Loop over all progenies and pairs of progenies */
- for (int j = 0; j < 8; j++) {
- if (ci->progeny[j] != NULL) {
- cell_activate_subcell_rt_tasks(ci->progeny[j], NULL, s);
- for (int k = j + 1; k < 8; k++)
- if (ci->progeny[k] != NULL)
- cell_activate_subcell_rt_tasks(ci->progeny[j], ci->progeny[k], s);
- }
- }
- } else {
- /* We have reached the bottom of the tree: activate drift */
- cell_activate_drift_spart(ci, s);
- cell_activate_drift_part(ci, s);
- }
- }
-
- /* Otherwise, pair interaction */
- else {
-
- /* Get the orientation of the pair. */
- double shift[3];
- const int sid = space_getsid(s->space, &ci, &cj, shift);
-
- const int ci_active =
- rt_should_iact_cell_pair(ci, cj, e) && (cj->hydro.count > 0);
- const int cj_active =
- rt_should_iact_cell_pair(cj, ci, e) && (ci->hydro.count > 0);
-
- /* Should we even bother? */
- if (!ci_active && !cj_active) return;
-
- /* recurse? */
- if (cell_can_recurse_in_pair_stars_task(ci, cj) &&
- cell_can_recurse_in_pair_stars_task(cj, ci)) {
-
- const struct cell_split_pair *csp = &cell_split_pairs[sid];
- for (int k = 0; k < csp->count; k++) {
- const int pid = csp->pairs[k].pid;
- const int pjd = csp->pairs[k].pjd;
- if (ci->progeny[pid] != NULL && cj->progeny[pjd] != NULL)
- cell_activate_subcell_rt_tasks(ci->progeny[pid], cj->progeny[pjd], s);
- }
- }
-
- /* Otherwise, activate the sorts and drifts. */
- else {
-
- if (ci_active) {
-
- /* We are going to interact this pair, so store some values. */
- atomic_or(&cj->hydro.requires_sorts, 1 << sid);
- atomic_or(&ci->stars.requires_sorts, 1 << sid);
-
- cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
- ci->stars.dx_max_sort_old = ci->stars.dx_max_sort;
-
- /* Activate the drifts if the cells are local. */
- if (ci->nodeID == engine_rank) cell_activate_drift_spart(ci, s);
- if (cj->nodeID == engine_rank) cell_activate_drift_part(cj, s);
-
- /* Do we need to sort the cells? */
- cell_activate_hydro_sorts(cj, sid, s);
- cell_activate_stars_sorts(ci, sid, s);
- }
-
- if (cj_active) {
-
- /* We are going to interact this pair, so store some values. */
- atomic_or(&cj->stars.requires_sorts, 1 << sid);
- atomic_or(&ci->hydro.requires_sorts, 1 << sid);
-
- ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
- cj->stars.dx_max_sort_old = cj->stars.dx_max_sort;
-
- /* Activate the drifts if the cells are local. */
- if (ci->nodeID == engine_rank) cell_activate_drift_part(ci, s);
- if (cj->nodeID == engine_rank) cell_activate_drift_spart(cj, s);
-
- /* Do we need to sort the cells? */
- cell_activate_hydro_sorts(ci, sid, s);
- cell_activate_stars_sorts(cj, sid, s);
- }
- }
- }
-}
-
/**
* @brief Un-skips all the hydro tasks associated with a given cell and checks
* if the space needs to be rebuilt.
@@ -2873,117 +2752,6 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s) {
const int nodeID = e->nodeID;
int rebuild = 0; /* TODO: implement rebuild conditions? */
- if (c->stars.drift != NULL) {
- if (rt_should_iact_cell(c, e)) {
- cell_activate_drift_part(c, s);
- cell_activate_drift_spart(c, s);
- }
- }
-
- /* Now unskip all RT specific interaction tasks */
- for (struct link *l = c->hydro.rt_inject; l != NULL; l = l->next) {
- struct task *t = l->t;
- struct cell *ci = t->ci;
- struct cell *cj = t->cj;
-#ifdef WITH_MPI
- const int ci_nodeID = ci->nodeID;
- const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
-#else
- const int ci_nodeID = nodeID;
- const int cj_nodeID = nodeID;
-#endif
-
- /* Activate the drifts */
- if (t->type == task_type_self) {
- if (rt_should_iact_cell(ci, e)) {
- cell_activate_drift_part(ci, s);
- cell_activate_drift_spart(ci, s);
- scheduler_activate(s, t);
- }
- }
-
- else if (t->type == task_type_pair) {
-
- const int ci_active = rt_should_iact_cell_pair(ci, cj, e);
- const int cj_active = (cj != NULL) && rt_should_iact_cell_pair(cj, ci, e);
-
- /* Only activate tasks that involve a local active cell. */
- if ((ci_active || cj_active) &&
- (ci_nodeID == nodeID || cj_nodeID == nodeID)) {
- scheduler_activate(s, t);
-
- /* Do ci */
- if (ci_active) {
- /* stars for ci */
- atomic_or(&ci->stars.requires_sorts, 1 << t->flags);
- ci->stars.dx_max_sort_old = ci->stars.dx_max_sort;
-
- /* hydro for cj */
- atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
- cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
-
- /* Activate the drift tasks. */
- if (ci_nodeID == nodeID) cell_activate_drift_spart(ci, s);
- if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
-
- /* Check the sorts and activate them if needed. */
- cell_activate_stars_sorts(ci, t->flags, s);
- cell_activate_hydro_sorts(cj, t->flags, s);
- }
-
- /* Do cj */
- if (cj_active) {
- /* hydro for ci */
- atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
- ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
-
- /* stars for cj */
- atomic_or(&cj->stars.requires_sorts, 1 << t->flags);
- cj->stars.dx_max_sort_old = cj->stars.dx_max_sort;
-
- /* Activate the drift tasks. */
- if (cj_nodeID == nodeID) cell_activate_drift_spart(cj, s);
- if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
-
- /* Check the sorts and activate them if needed. */
- cell_activate_hydro_sorts(ci, t->flags, s);
- cell_activate_stars_sorts(cj, t->flags, s);
- }
- }
- }
-
- else if (t->type == task_type_sub_self) {
- scheduler_activate(s, t);
- cell_activate_subcell_rt_tasks(ci, NULL, s);
- }
-
- else if (t->type == task_type_sub_pair) {
- scheduler_activate(s, t);
- cell_activate_subcell_rt_tasks(ci, cj, s);
- }
-
- /* Only interested in pair interactions as of here. */
- if (t->type == task_type_pair || t->type == task_type_sub_pair) {
- /* Check whether there was too much particle motion, i.e. the */
- /* cell neighbour conditions were violated. */
- if (cell_need_rebuild_for_stars_pair(ci, cj)) rebuild = 1;
- if (cell_need_rebuild_for_stars_pair(cj, ci)) rebuild = 1;
-
- /* Activate rt_in for each cell that is part of
- * a pair/sub_pair task as to not miss any dependencies */
- if (ci_nodeID == nodeID)
- scheduler_activate(s, ci->hydro.super->hydro.rt_in);
- if (cj_nodeID == nodeID)
- scheduler_activate(s, cj->hydro.super->hydro.rt_in);
-
- /* For the same reason, catch the dependencies with the RT ghost1 */
- if (ci_nodeID == nodeID)
- scheduler_activate(s, ci->hydro.super->hydro.rt_ghost1);
- if (cj_nodeID == nodeID)
- scheduler_activate(s, cj->hydro.super->hydro.rt_ghost1);
- }
- }
-
if (c->nodeID == nodeID) {
for (struct link *l = c->hydro.rt_gradient; l != NULL; l = l->next) {
@@ -3003,7 +2771,7 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s) {
#endif
const int ci_active = cell_is_active_hydro(ci, e);
- const int cj_active = (cj != NULL) && cell_is_active_hydro(cj, e);
+ const int cj_active = ((cj != NULL) && cell_is_active_hydro(cj, e));
if (t->type == task_type_self || t->type == task_type_sub_self) {
if (ci_active) scheduler_activate(s, t);
@@ -3043,7 +2811,7 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s) {
#endif
const int ci_active = cell_is_active_hydro(ci, e);
- const int cj_active = (cj != NULL) && cell_is_active_hydro(cj, e);
+ const int cj_active = ((cj != NULL) && cell_is_active_hydro(cj, e));
if (t->type == task_type_self || t->type == task_type_sub_self) {
if (ci_active) scheduler_activate(s, t);
@@ -3068,22 +2836,9 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s) {
/* Unskip all the other task types */
- if (rt_should_do_unskip_cell(c, e)) {
- /* If we don't have any pair/sub_pair tasks for this cell, then we haven't
- * unskipped the ghosts yet. Do this now. */
-
- /* You need to pay attention to stars as well when unskipping rt_in
- * to gather dependencies from the feedback loop. */
+ if (cell_is_active_hydro(c, e)) {
if (c->hydro.rt_in != NULL) scheduler_activate(s, c->hydro.rt_in);
-
- /* Also activate the rt_ghost1 task even if you don't have stars in
- * this cell to gather dependencies properly. Otherwise, dependency
- * issues arise when the timestep task starts changing what's active
- * and what's not active */
if (c->hydro.rt_ghost1 != NULL) scheduler_activate(s, c->hydro.rt_ghost1);
- }
-
- if (cell_is_active_hydro(c, e)) {
if (c->hydro.rt_ghost2 != NULL) scheduler_activate(s, c->hydro.rt_ghost2);
if (c->hydro.rt_transport_out != NULL)
scheduler_activate(s, c->hydro.rt_transport_out);
diff --git a/src/engine.c b/src/engine.c
index 7ad7323ebf1c586a974f42e36bfeca8bc37103e5..6b9e7b094fa21a302b0ebaa725dc93a4a982f5c4 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1119,12 +1119,11 @@ int engine_estimate_nr_tasks(const struct engine *e) {
}
#endif
if (e->policy & engine_policy_rt) {
- /* inject: 1 self + (3^3-1)/2 = 26/2 = 13 pairs | 14
- * gradient: 1 self + 13 pairs | + 14
+ /* gradient: 1 self + 13 pairs | 14
* transport: 1 self + 13 pairs | + 14
* implicits: in + out, transport_out | + 3
* others: ghost1, ghost2, thermochemistry | + 3 */
- n1 += 48;
+ n1 += 34;
}
#ifdef WITH_MPI
@@ -1627,7 +1626,6 @@ void engine_skip_force_and_kick(struct engine *e) {
t->subtype == task_subtype_part_prep1 ||
t->subtype == task_subtype_spart_prep2 ||
t->subtype == task_subtype_sf_counts ||
- t->subtype == task_subtype_rt_inject ||
t->subtype == task_subtype_rt_gradient ||
t->subtype == task_subtype_rt_transport)
t->skip = 1;
@@ -1929,8 +1927,10 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
TIMER_TOC2(timer_runners);
/* Initialise additional RT data now that time bins are set */
+#ifdef SWIFT_RT_DEBUG_CHECKS
if (e->policy & engine_policy_rt)
space_convert_rt_quantities_after_zeroth_step(e->s, e->verbose);
+#endif
#ifdef SWIFT_HYDRO_DENSITY_CHECKS
/* Run the brute-force hydro calculation for some parts */
@@ -2439,12 +2439,6 @@ void engine_step(struct engine *e) {
space_check_unskip_flags(e->s);
#endif
-#if defined(SWIFT_RT_DEBUG_CHECKS)
- /* if we're running the debug RT scheme, do some checks after every step */
- if (e->policy & engine_policy_rt)
- rt_debugging_checks_end_of_step(e, e->verbose);
-#endif
-
/* Compute the local accumulated deadtime. */
const ticks deadticks = (e->nr_threads * e->sched.deadtime.waiting_ticks) -
e->sched.deadtime.active_ticks;
@@ -2482,6 +2476,13 @@ void engine_step(struct engine *e) {
engine_check_for_dumps(e);
+#if defined(SWIFT_RT_DEBUG_CHECKS)
+ /* if we're running the debug RT scheme, do some checks after every step.
+ * Do this after the output so we can safely reset debugging checks now. */
+ if (e->policy & engine_policy_rt)
+ rt_debugging_checks_end_of_step(e, e->verbose);
+#endif
+
TIMER_TOC2(timer_step);
clocks_gettime(&time2);
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index ce9158d155398d9dc74d7077737b3064a05c3a40..3397544b764b7a6b6aa7000ccd8d0bb3ce77e996 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -1491,12 +1491,7 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
/* non-implicit ghost 1 */
c->hydro.rt_ghost1 = scheduler_addtask(
s, task_type_rt_ghost1, task_subtype_none, 0, 0, c, NULL);
- /* add the explicit dependency on kick2 for cases where injection
- * gets skipped */
- scheduler_addunlock(s, c->super->kick2, c->hydro.rt_ghost1);
- /* add the explicit dependency for the timestep task for the
- * cases where we have active sparts, but no active parts */
- scheduler_addunlock(s, c->hydro.rt_ghost1, c->super->timestep);
+ scheduler_addunlock(s, c->hydro.rt_in, c->hydro.rt_ghost1);
/* non-implicit ghost 2 */
c->hydro.rt_ghost2 = scheduler_addtask(
@@ -2155,7 +2150,6 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
struct task *t_do_gas_swallow = NULL;
struct task *t_do_bh_swallow = NULL;
struct task *t_bh_feedback = NULL;
- struct task *t_rt_inject = NULL;
struct task *t_sink_formation = NULL;
struct task *t_rt_gradient = NULL;
struct task *t_rt_transport = NULL;
@@ -2256,8 +2250,6 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
if (with_rt) {
- t_rt_inject = scheduler_addtask(
- sched, task_type_self, task_subtype_rt_inject, flags, 0, ci, NULL);
t_rt_gradient =
scheduler_addtask(sched, task_type_self, task_subtype_rt_gradient,
flags, 0, ci, NULL);
@@ -2292,7 +2284,6 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
}
if (with_rt) {
- engine_addlink(e, &ci->hydro.rt_inject, t_rt_inject);
engine_addlink(e, &ci->hydro.rt_gradient, t_rt_gradient);
engine_addlink(e, &ci->hydro.rt_transport, t_rt_transport);
}
@@ -2430,11 +2421,6 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
if (with_rt) {
- scheduler_addunlock(sched, ci->hydro.super->stars.drift, t_rt_inject);
- scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t_rt_inject);
- scheduler_addunlock(sched, ci->hydro.super->hydro.rt_in, t_rt_inject);
- scheduler_addunlock(sched, t_rt_inject,
- ci->hydro.super->hydro.rt_ghost1);
scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t_rt_gradient);
scheduler_addunlock(sched, ci->hydro.super->hydro.rt_ghost1,
t_rt_gradient);
@@ -2522,8 +2508,6 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
if (with_rt) {
- t_rt_inject = scheduler_addtask(
- sched, task_type_pair, task_subtype_rt_inject, flags, 0, ci, cj);
t_rt_gradient = scheduler_addtask(
sched, task_type_pair, task_subtype_rt_gradient, flags, 0, ci, cj);
t_rt_transport = scheduler_addtask(
@@ -2572,8 +2556,6 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
engine_addlink(e, &cj->black_holes.feedback, t_bh_feedback);
}
if (with_rt) {
- engine_addlink(e, &ci->hydro.rt_inject, t_rt_inject);
- engine_addlink(e, &cj->hydro.rt_inject, t_rt_inject);
engine_addlink(e, &ci->hydro.rt_gradient, t_rt_gradient);
engine_addlink(e, &cj->hydro.rt_gradient, t_rt_gradient);
engine_addlink(e, &ci->hydro.rt_transport, t_rt_transport);
@@ -2637,11 +2619,9 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
}
if (with_rt) {
- scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t_rt_inject);
scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t_rt_gradient);
if (ci->hydro.super != cj->hydro.super) {
- scheduler_addunlock(sched, cj->hydro.super->hydro.sorts, t_rt_inject);
scheduler_addunlock(sched, cj->hydro.super->hydro.sorts,
t_rt_gradient);
}
@@ -2768,12 +2748,6 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
if (with_rt) {
- scheduler_addunlock(sched, ci->hydro.super->stars.sorts, t_rt_inject);
- scheduler_addunlock(sched, ci->hydro.super->stars.drift, t_rt_inject);
- scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t_rt_inject);
- scheduler_addunlock(sched, ci->hydro.super->hydro.rt_in, t_rt_inject);
- scheduler_addunlock(sched, t_rt_inject,
- ci->hydro.super->hydro.rt_ghost1);
scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
t_rt_gradient);
scheduler_addunlock(sched, ci->hydro.super->hydro.rt_ghost1,
@@ -2914,16 +2888,6 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
if (with_rt) {
- scheduler_addunlock(sched, cj->hydro.super->stars.sorts,
- t_rt_inject);
- scheduler_addunlock(sched, cj->hydro.super->stars.drift,
- t_rt_inject);
- scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
- t_rt_inject);
- scheduler_addunlock(sched, cj->hydro.super->hydro.rt_in,
- t_rt_inject);
- scheduler_addunlock(sched, t_rt_inject,
- cj->hydro.super->hydro.rt_ghost1);
scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
t_rt_gradient);
scheduler_addunlock(sched, cj->hydro.super->hydro.rt_ghost1,
@@ -3051,9 +3015,6 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
if (with_rt) {
- t_rt_inject =
- scheduler_addtask(sched, task_type_sub_self, task_subtype_rt_inject,
- flags, 0, ci, NULL);
t_rt_gradient =
scheduler_addtask(sched, task_type_sub_self,
task_subtype_rt_gradient, flags, 0, ci, NULL);
@@ -3088,7 +3049,6 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
}
if (with_rt) {
- engine_addlink(e, &ci->hydro.rt_inject, t_rt_inject);
engine_addlink(e, &ci->hydro.rt_gradient, t_rt_gradient);
engine_addlink(e, &ci->hydro.rt_transport, t_rt_transport);
}
@@ -3235,13 +3195,6 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
if (with_rt) {
- scheduler_addunlock(sched, ci->hydro.super->stars.drift, t_rt_inject);
- scheduler_addunlock(sched, ci->hydro.super->stars.sorts, t_rt_inject);
- scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t_rt_inject);
- scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t_rt_inject);
- scheduler_addunlock(sched, ci->hydro.super->hydro.rt_in, t_rt_inject);
- scheduler_addunlock(sched, t_rt_inject,
- ci->hydro.super->hydro.rt_ghost1);
scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t_rt_gradient);
scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t_rt_gradient);
scheduler_addunlock(sched, ci->hydro.super->hydro.rt_ghost1,
@@ -3337,9 +3290,6 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
if (with_rt) {
- t_rt_inject =
- scheduler_addtask(sched, task_type_sub_pair, task_subtype_rt_inject,
- flags, 0, ci, cj);
t_rt_gradient =
scheduler_addtask(sched, task_type_sub_pair,
task_subtype_rt_gradient, flags, 0, ci, cj);
@@ -3392,8 +3342,6 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
engine_addlink(e, &cj->black_holes.feedback, t_bh_feedback);
}
if (with_rt) {
- engine_addlink(e, &ci->hydro.rt_inject, t_rt_inject);
- engine_addlink(e, &cj->hydro.rt_inject, t_rt_inject);
engine_addlink(e, &ci->hydro.rt_gradient, t_rt_gradient);
engine_addlink(e, &cj->hydro.rt_gradient, t_rt_gradient);
engine_addlink(e, &ci->hydro.rt_transport, t_rt_transport);
@@ -3457,10 +3405,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
if (with_rt) {
- scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t_rt_inject);
scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t_rt_gradient);
if (ci->hydro.super != cj->hydro.super) {
- scheduler_addunlock(sched, cj->hydro.super->hydro.sorts, t_rt_inject);
scheduler_addunlock(sched, cj->hydro.super->hydro.sorts,
t_rt_gradient);
}
@@ -3587,12 +3533,6 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
if (with_rt) {
- scheduler_addunlock(sched, ci->hydro.super->stars.sorts, t_rt_inject);
- scheduler_addunlock(sched, ci->hydro.super->stars.drift, t_rt_inject);
- scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t_rt_inject);
- scheduler_addunlock(sched, ci->hydro.super->hydro.rt_in, t_rt_inject);
- scheduler_addunlock(sched, t_rt_inject,
- ci->hydro.super->hydro.rt_ghost1);
scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
t_rt_gradient);
scheduler_addunlock(sched, ci->hydro.super->hydro.rt_ghost1,
@@ -3731,16 +3671,6 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
cj->hydro.super->black_holes.black_holes_out);
}
if (with_rt) {
- scheduler_addunlock(sched, cj->hydro.super->stars.sorts,
- t_rt_inject);
- scheduler_addunlock(sched, cj->hydro.super->stars.drift,
- t_rt_inject);
- scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
- t_rt_inject);
- scheduler_addunlock(sched, cj->hydro.super->hydro.rt_in,
- t_rt_inject);
- scheduler_addunlock(sched, t_rt_inject,
- cj->hydro.super->hydro.rt_ghost1);
scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
t_rt_gradient);
scheduler_addunlock(sched, cj->hydro.super->hydro.rt_ghost1,
diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c
index c1cc4567209b640ae5def59aba769ae174652b6f..d01e6981c60e74ccdad14620ea7ba5a4cb2a7a54 100644
--- a/src/engine_marktasks.c
+++ b/src/engine_marktasks.c
@@ -52,7 +52,6 @@
#include "error.h"
#include "feedback.h"
#include "proxy.h"
-#include "rt_active.h"
#include "timers.h"
/**
@@ -103,7 +102,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
cell_is_active_sinks(ci, e) || ci_active_hydro;
const int ci_active_stars = cell_need_activating_stars(
ci, e, with_star_formation, with_star_formation_sink);
- const int ci_active_rt = with_rt && rt_should_do_unskip_cell(ci, e);
+ const int ci_active_rt = ci_active_hydro && with_rt;
/* Activate the hydro drift */
if (t_type == task_type_self && t_subtype == task_subtype_density) {
@@ -347,26 +346,9 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
/* Activate RT tasks */
- else if (t_type == task_type_self &&
- t_subtype == task_subtype_rt_inject) {
- if (ci_active_rt) {
- scheduler_activate(s, t);
- cell_activate_drift_part(ci, s);
- cell_activate_drift_spart(ci, s);
- }
- }
-
- else if (t_type == task_type_sub_self &&
- t_subtype == task_subtype_rt_inject) {
- if (ci_active_rt) {
- scheduler_activate(s, t);
- cell_activate_subcell_rt_tasks(ci, NULL, s);
- }
- }
-
else if (t_subtype == task_subtype_rt_gradient ||
t_subtype == task_subtype_rt_transport) {
- if (ci_active_hydro) scheduler_activate(s, t);
+ if (ci_active_rt) scheduler_activate(s, t);
}
#ifdef SWIFT_DEBUG_CHECKS
@@ -408,8 +390,8 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
const int cj_active_stars = cell_need_activating_stars(
cj, e, with_star_formation, with_star_formation_sink);
- const int ci_active_rt = with_rt && rt_should_iact_cell_pair(ci, cj, e);
- const int cj_active_rt = with_rt && rt_should_iact_cell_pair(cj, ci, e);
+ const int ci_active_rt = ci_active_hydro && with_rt;
+ const int cj_active_rt = cj_active_hydro && with_rt;
/* Only activate tasks that involve a local active cell. */
if ((t_subtype == task_subtype_density ||
@@ -780,99 +762,12 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
}
- /* RT injection tasks */
- else if (t_subtype == task_subtype_rt_inject) {
-
- /* We only want to activate the task if the cell is active and is
- going to update some gas on the *local* node */
- if ((ci_nodeID == nodeID && cj_nodeID == nodeID) &&
- (ci_active_rt || cj_active_rt)) {
-
- scheduler_activate(s, t);
-
- /* Set the correct sorting flags */
- if (t_type == task_type_pair) {
-
- /* Add rt_in dependencies for each cell that is part of
- * a pair task as to not miss any dependencies */
- if (ci_nodeID == nodeID)
- scheduler_activate(s, ci->hydro.super->hydro.rt_in);
- if (cj_nodeID == nodeID)
- scheduler_activate(s, cj->hydro.super->hydro.rt_in);
-
- /* For the same reason, catch the dependencies with the RT ghost1 */
- if (ci_nodeID == nodeID)
- scheduler_activate(s, ci->hydro.super->hydro.rt_ghost1);
- if (cj_nodeID == nodeID)
- scheduler_activate(s, cj->hydro.super->hydro.rt_ghost1);
-
- /* Do ci */
- if (ci_active_rt) {
-
- /* stars for ci */
- atomic_or(&ci->stars.requires_sorts, 1 << t->flags);
- ci->stars.dx_max_sort_old = ci->stars.dx_max_sort;
-
- /* hydro for cj */
- atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
- cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
-
- /* Activate the drift tasks. */
- if (ci_nodeID == nodeID) cell_activate_drift_spart(ci, s);
- if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
-
- /* Check the sorts and activate them if needed. */
- cell_activate_hydro_sorts(cj, t->flags, s);
- cell_activate_stars_sorts(ci, t->flags, s);
- }
-
- /* Do cj */
- if (cj_active_rt) {
-
- /* hydro for ci */
- atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
- ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
-
- /* stars for cj */
- atomic_or(&cj->stars.requires_sorts, 1 << t->flags);
- cj->stars.dx_max_sort_old = cj->stars.dx_max_sort;
-
- /* Activate the drift tasks. */
- if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
- if (cj_nodeID == nodeID) cell_activate_drift_spart(cj, s);
-
- /* Check the sorts and activate them if needed. */
- cell_activate_hydro_sorts(ci, t->flags, s);
- cell_activate_stars_sorts(cj, t->flags, s);
- }
- }
-
- /* Store current values of dx_max and h_max. */
- else if (t_type == task_type_sub_pair) {
- cell_activate_subcell_rt_tasks(ci, cj, s);
-
- /* Add rt_in dependencies for each cell that is part of
- * a sub_pair task as to not miss any dependencies */
- if (ci_nodeID == nodeID)
- scheduler_activate(s, ci->hydro.super->hydro.rt_in);
- if (cj_nodeID == nodeID)
- scheduler_activate(s, cj->hydro.super->hydro.rt_in);
-
- /* For the same reason, catch the dependencies with the RT ghost1 */
- if (ci_nodeID == nodeID)
- scheduler_activate(s, ci->hydro.super->hydro.rt_ghost1);
- if (cj_nodeID == nodeID)
- scheduler_activate(s, cj->hydro.super->hydro.rt_ghost1);
- }
- }
- }
-
/* RT gradient and transport tasks */
else if (t_subtype == task_subtype_rt_gradient) {
/* We only want to activate the task if the cell is active and is
going to update some gas on the *local* node */
if ((ci_nodeID == nodeID && cj_nodeID == nodeID) &&
- (ci_active_hydro || cj_active_hydro)) {
+ (ci_active_rt || cj_active_rt)) {
/* The gradient and transport task subtypes mirror the hydro tasks.
* Therefore all the (subcell) sorts and drifts should already have
* been activated properly in the hydro part of the activation. */
@@ -893,7 +788,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
/* We only want to activate the task if the cell is active and is
going to update some gas on the *local* node */
if ((ci_nodeID == nodeID && cj_nodeID == nodeID) &&
- (ci_active_hydro || cj_active_hydro)) {
+ (ci_active_rt || cj_active_rt)) {
/* The gradient and transport task subtypes mirror the hydro tasks.
* Therefore all the (subcell) sorts and drifts should already have
* been activated properly in the hydro part of the activation. */
@@ -1284,8 +1179,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
#endif
} /* Only interested in RT tasks as of here. */
#ifdef WITH_MPI
- else if (t_subtype == task_subtype_rt_inject ||
- t_subtype == task_subtype_rt_gradient ||
+ else if (t_subtype == task_subtype_rt_gradient ||
t_subtype == task_subtype_rt_transport) {
error("RT doesn't work with MPI yet.");
}
@@ -1443,22 +1337,15 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
/* Radiative transfer implicit tasks */
else if (t->type == task_type_rt_in) {
- if (rt_should_do_unskip_cell(t->ci, e)) scheduler_activate(s, t);
- }
-
- else if (t->type == task_type_rt_transport_out ||
- t->type == task_type_rt_out) {
- if (cell_is_active_hydro(t->ci, e)) scheduler_activate(s, t);
- }
-
- /* rt_ghost1 is special: Also check for stars activity to catch
- * dependencies further down the line (e.g. timestep task) */
- else if (t->type == task_type_rt_ghost1) {
- if (rt_should_do_unskip_cell(t->ci, e)) scheduler_activate(s, t);
+ if (cell_is_active_hydro(t->ci, e) ||
+ cell_need_activating_stars(t->ci, e, with_star_formation,
+ with_star_formation_sink))
+ scheduler_activate(s, t);
}
- /* Radiative transfer ghosts and thermochemistry*/
- else if (t->type == task_type_rt_ghost2 || t->type == task_type_rt_tchem) {
+ else if (t->type == task_type_rt_ghost1 || t->type == task_type_rt_ghost2 ||
+ t->type == task_type_rt_transport_out ||
+ t->type == task_type_rt_tchem || t->type == task_type_rt_out) {
if (cell_is_active_hydro(t->ci, e)) scheduler_activate(s, t);
}
diff --git a/src/engine_unskip.c b/src/engine_unskip.c
index eb81c2ad42fb03b970054b2b4112fde79e4f907c..78c18a3e470b7670b9d2d438c328679c346ed3bd 100644
--- a/src/engine_unskip.c
+++ b/src/engine_unskip.c
@@ -28,7 +28,6 @@
#include "active.h"
#include "cell.h"
#include "memswap.h"
-#include "rt_active.h"
/* Load the profiler header, if needed. */
#ifdef WITH_PROFILER
@@ -261,7 +260,7 @@ static void engine_do_unskip_rt(struct cell *c, struct engine *e) {
if (!cell_get_flag(c, cell_flag_has_tasks)) return;
/* Do we have work to do? */
- if (!rt_should_do_unskip_cell(c, e)) return;
+ if (!cell_is_active_hydro(c, e)) return;
/* Recurse */
if (c->split) {
@@ -417,7 +416,7 @@ void engine_unskip(struct engine *e) {
(with_stars && c->nodeID == nodeID && cell_is_active_stars(c, e)) ||
(with_sinks && cell_is_active_sinks(c, e)) ||
(with_black_holes && cell_is_active_black_holes(c, e)) ||
- (with_rt && rt_should_do_unskip_cell(c, e))) {
+ (with_rt && cell_is_active_hydro(c, e))) {
if (num_active_cells != k)
memswap(&local_cells[k], &local_cells[num_active_cells], sizeof(int));
diff --git a/src/rt/GEAR/rt.h b/src/rt/GEAR/rt.h
index a53a8287de7acc93825fe656a7a5377466c36452..842782719394e25522f7e64e8ba9c5e74344ae95 100644
--- a/src/rt/GEAR/rt.h
+++ b/src/rt/GEAR/rt.h
@@ -56,6 +56,10 @@ rt_compute_stellar_emission_rate(struct spart* restrict sp, double time,
const struct phys_const* phys_const,
const struct unit_system* internal_units) {
+#ifdef SWIFT_RT_DEBUG_CHECKS
+ sp->rt_data.debug_emission_rate_set += 1;
+#endif
+
/* Skip initial fake time-step */
if (dt == 0.0l) return;
@@ -98,10 +102,7 @@ __attribute__((always_inline)) INLINE static void rt_reset_part(
/* reset this here as well as in the rt_debugging_checks_end_of_step()
* routine to test task dependencies are done right */
p->rt_data.debug_iact_stars_inject = 0;
- p->rt_data.debug_iact_stars_inject_prep = 0;
- /* skip this for GEAR */
- /* p->rt_data.debug_injection_check = 0; */
p->rt_data.debug_calls_iact_gradient_interaction = 0;
p->rt_data.debug_calls_iact_transport_interaction = 0;
@@ -134,14 +135,13 @@ __attribute__((always_inline)) INLINE static void rt_first_init_part(
#ifdef SWIFT_RT_DEBUG_CHECKS
p->rt_data.debug_radiation_absorbed_tot = 0ULL;
- p->rt_data.debug_iact_stars_inject_prep_tot = 0ULL;
#endif
}
/**
* @brief Initialises particle quantities that can't be set
* otherwise before the zeroth step is finished. E.g. because
- * they require the particle density to be known.
+ * they require the particle density and time step to be known.
*
* @param p particle to work on
* @param rt_props RT properties struct
@@ -154,7 +154,14 @@ rt_init_part_after_zeroth_step(struct part* restrict p,
/* If we're running with debugging checks on, reset debugging
* counters and flags in particular after the zeroth step so
* that the checks work as intended. */
+ rt_init_part(p);
rt_reset_part(p);
+ /* Since the inject_prep has been moved to the density loop, the
+ * initialization at startup is messing with the total counters for stars
+ * because the density is called, but not the force-and-kick tasks. So reset
+ * the total counters here as well so that they will match the star counters.
+ */
+ p->rt_data.debug_radiation_absorbed_tot = 0ULL;
#endif
}
@@ -171,6 +178,22 @@ __attribute__((always_inline)) INLINE static void rt_init_spart(
for (int i = 0; i < 8; i++) {
sp->rt_data.octant_weights[i] = 0.f;
}
+
+#ifdef SWIFT_RT_DEBUG_CHECKS
+ /* reset this here as well as in the rt_debugging_checks_end_of_step()
+ * routine to test task dependencies are done right */
+ sp->rt_data.debug_iact_hydro_inject_prep = 0;
+ sp->rt_data.debug_iact_hydro_inject = 0;
+ sp->rt_data.debug_emission_rate_set = 0;
+
+ for (int g = 0; g < RT_NGROUPS; g++) {
+ sp->rt_data.debug_injected_energy[g] = 0.f;
+ }
+ for (int g = 0; g < RT_NGROUPS; g++) {
+ sp->rt_data.emission_this_step[g] = 0.f;
+ }
+ sp->rt_data.debug_psi_sum = 0.f;
+#endif
}
/**
@@ -188,21 +211,6 @@ __attribute__((always_inline)) INLINE static void rt_reset_spart(
for (int g = 0; g < RT_NGROUPS; g++) {
sp->rt_data.emission_this_step[g] = 0.f;
}
-
-#ifdef SWIFT_RT_DEBUG_CHECKS
- /* reset this here as well as in the rt_debugging_checks_end_of_step()
- * routine to test task dependencies are done right */
- sp->rt_data.debug_iact_hydro_inject = 0;
- sp->rt_data.debug_iact_hydro_inject_prep = 0;
-
- sp->rt_data.debug_emission_rate_set = 0;
- /* skip this for GEAR */
- /* sp->rt_data.debug_injection_check = 0; */
-
- for (int g = 0; g < RT_NGROUPS; g++) {
- sp->rt_data.debug_injected_energy[g] = 0.f;
- }
-#endif
}
/**
@@ -226,7 +234,7 @@ __attribute__((always_inline)) INLINE static void rt_first_init_spart(
/**
* @brief Initialises particle quantities that can't be set
* otherwise before the zeroth step is finished. E.g. because
- * they require the star density to be known.
+ * they require the star density and time step to be known.
* @param sp star particle to work on
* @param time current system time
* @param star_age age of the star *at the end of the step*
@@ -246,14 +254,14 @@ rt_init_star_after_zeroth_step(struct spart* restrict sp, double time,
/* If we're running with debugging checks on, reset debugging
* counters and flags in particular after the zeroth step so
* that the checks work as intended. */
+ rt_init_spart(sp);
rt_reset_spart(sp);
+ /* Since the inject_prep has been moved to the density loop, the
+ * initialization at startup is messing with the total counters because
+ * the density is called, but not the force-and-kick tasks. So reset
+ * the total counters here as well. */
+ sp->rt_data.debug_radiation_emitted_tot = 0ULL;
#endif
-
- /* If we're running with star controlled injection, we don't
- * need to compute the stellar emission rates here now. */
- if (rt_props->hydro_controlled_injection)
- rt_compute_stellar_emission_rate(sp, time, star_age, dt, rt_props,
- phys_const, internal_units);
}
/**
@@ -280,7 +288,7 @@ __attribute__((always_inline)) INLINE static void rt_part_has_no_neighbours(
/**
* @brief Exception handle a star part not having any neighbours in ghost task
*
- * @param sp The star particle to work on
+ * @param sp The #spart.
*/
__attribute__((always_inline)) INLINE static void rt_spart_has_no_neighbours(
struct spart* sp) {
@@ -530,6 +538,20 @@ __attribute__((always_inline)) INLINE static void rt_tchem(
const struct phys_const* restrict phys_const,
const struct unit_system* restrict us, const double dt) {
+#ifdef SWIFT_RT_DEBUG_CHECKS
+ if (p->rt_data.debug_kicked != 1)
+ error("Trying to do thermochemistry on unkicked particle %lld (count=%d)",
+ p->id, p->rt_data.debug_kicked);
+ if (!p->rt_data.debug_injection_done)
+ error("Trying to do thermochemistry when injection step hasn't been done");
+ if (!p->rt_data.debug_gradients_done)
+ error("Trying to do thermochemistry when gradient step hasn't been done");
+ if (!p->rt_data.debug_transport_done)
+ error("Trying to do thermochemistry when transport step hasn't been done");
+
+ p->rt_data.debug_thermochem_done += 1;
+#endif
+
/* Note: Can't pass rt_props as const struct because of grackle
* accessinging its properties there */
@@ -632,20 +654,25 @@ __attribute__((always_inline)) INLINE static void rt_prepare_force(
* @brief Clean the allocated memory inside the RT properties struct.
*
* @param props the #rt_props.
+ * @param restart did we restart?
*/
__attribute__((always_inline)) INLINE static void rt_clean(
- struct rt_props* props) {
-
- /* Clean up grackle data. This is a call to a grackle function */
- _free_chemistry_data(&props->grackle_chemistry_data,
- props->grackle_chemistry_rates);
-
- for (int g = 0; g < RT_NGROUPS; g++) {
- free(props->energy_weighted_cross_sections[g]);
- free(props->number_weighted_cross_sections[g]);
+ struct rt_props* props, int restart) {
+
+ /* If we were restarting, free-ing manually will lead to
+ * segfaults since we didn't malloc the stuff */
+ if (!restart) {
+ /* Clean up grackle data. This is a call to a grackle function */
+ _free_chemistry_data(&props->grackle_chemistry_data,
+ props->grackle_chemistry_rates);
+
+ for (int g = 0; g < RT_NGROUPS; g++) {
+ free(props->energy_weighted_cross_sections[g]);
+ free(props->number_weighted_cross_sections[g]);
+ }
+ free(props->energy_weighted_cross_sections);
+ free(props->number_weighted_cross_sections);
}
- free(props->energy_weighted_cross_sections);
- free(props->number_weighted_cross_sections);
#ifdef SWIFT_RT_DEBUG_CHECKS
fclose(props->conserved_energy_filep);
diff --git a/src/rt/GEAR/rt_debugging.h b/src/rt/GEAR/rt_debugging.h
index c6bd4d163fc960d9b905dc60fc1f71598425d5ea..5421dde78016efa4c1a512108d6389429bf46dd6 100644
--- a/src/rt/GEAR/rt_debugging.h
+++ b/src/rt/GEAR/rt_debugging.h
@@ -39,24 +39,18 @@ static void rt_debugging_end_of_step_stars_mapper(void *restrict map_data,
struct spart *restrict sparts = (struct spart *)map_data;
const struct engine *restrict e = (struct engine *)extra_data;
- int emission_sum = 0;
- int iacts_with_parts_sum = 0;
+ unsigned long long emission_sum_this_step = 0ULL;
unsigned long long emission_sum_tot = 0ULL;
- unsigned long long iacts_with_parts_sum_tot = 0ULL;
float emitted_energy[RT_NGROUPS];
for (int g = 0; g < RT_NGROUPS; g++) emitted_energy[g] = 0.f;
for (int k = 0; k < scount; k++) {
struct spart *restrict sp = &sparts[k];
- emission_sum += sp->rt_data.debug_iact_hydro_inject;
+ emission_sum_this_step += sp->rt_data.debug_iact_hydro_inject;
emission_sum_tot += sp->rt_data.debug_radiation_emitted_tot;
/* Reset all values here in case stars won't be active next step */
sp->rt_data.debug_iact_hydro_inject = 0;
-
- iacts_with_parts_sum += sp->rt_data.debug_iact_hydro_inject_prep;
- iacts_with_parts_sum_tot += sp->rt_data.debug_iact_hydro_inject_prep_tot;
- /* Reset all values here in case stars won't be active next step */
sp->rt_data.debug_iact_hydro_inject_prep = 0;
for (int g = 0; g < RT_NGROUPS; g++) {
@@ -70,12 +64,22 @@ static void rt_debugging_end_of_step_stars_mapper(void *restrict map_data,
float diff = 1.f - sp->rt_data.emission_this_step[g] /
sp->rt_data.debug_injected_energy[g];
- if (fabs(diff) > 1e-3) {
- message(
- "Incorrect injection ID %lld: "
- "group %d expected %.3g got %.3g diff %.3g",
- sp->id, g, sp->rt_data.emission_this_step[g],
- sp->rt_data.debug_injected_energy[g], diff);
+ if (fabsf(diff) > 1e-4) {
+ /* Dividing the total into several parts and summing them up again
+ * while hoping to obtain the same results may lead to diappointment
+ * due to roundoff errors. Check that the sum of the individual
+ * weights and the ones we collected for the injection are close
+ * enough. */
+ float psi_sum_now = 0.f;
+ for (int i = 0; i < 8; i++)
+ psi_sum_now += sp->rt_data.octant_weights[i];
+ float diff_weights = 1.f - sp->rt_data.debug_psi_sum / psi_sum_now;
+ if (fabsf(diff_weights) > 1e-4)
+ message(
+ "Incorrect injection ID %lld: "
+ "group %d expected %.3g got %.3g diff %.3g diff_weights %.3g",
+ sp->id, g, sp->rt_data.emission_this_step[g],
+ sp->rt_data.debug_injected_energy[g], diff, diff_weights);
}
}
emitted_energy[g] += sp->rt_data.debug_injected_energy[g];
@@ -84,14 +88,14 @@ static void rt_debugging_end_of_step_stars_mapper(void *restrict map_data,
for (int g = 0; g < RT_NGROUPS; g++) {
sp->rt_data.debug_injected_energy[g] = 0.f;
}
+ for (int g = 0; g < RT_NGROUPS; g++) {
+ sp->rt_data.emission_this_step[g] = 0.f;
+ }
}
- atomic_add(&e->rt_props->debug_radiation_emitted_this_step, emission_sum);
+ atomic_add(&e->rt_props->debug_radiation_emitted_this_step,
+ emission_sum_this_step);
atomic_add(&e->rt_props->debug_radiation_emitted_tot, emission_sum_tot);
- atomic_add(&e->rt_props->debug_star_injection_prep_iacts_with_parts_this_step,
- iacts_with_parts_sum);
- atomic_add(&e->rt_props->debug_star_injection_prep_iacts_with_parts_tot,
- iacts_with_parts_sum_tot);
for (int g = 0; g < RT_NGROUPS; g++)
atomic_add_f(&e->rt_props->debug_total_star_emitted_energy[g],
emitted_energy[g]);
@@ -107,26 +111,19 @@ static void rt_debugging_end_of_step_hydro_mapper(void *restrict map_data,
struct part *restrict parts = (struct part *)map_data;
const struct engine *restrict e = (struct engine *)extra_data;
- int absorption_sum = 0;
- int iacts_with_stars_sum = 0;
+ unsigned long long absorption_sum_this_step = 0ULL;
unsigned long long absorption_sum_tot = 0ULL;
- unsigned long long iacts_with_stars_sum_tot = 0ULL;
float energy_sum[RT_NGROUPS];
for (int g = 0; g < RT_NGROUPS; g++) energy_sum[g] = 0.f;
for (int k = 0; k < count; k++) {
struct part *restrict p = &parts[k];
- absorption_sum += p->rt_data.debug_iact_stars_inject;
+ absorption_sum_this_step += p->rt_data.debug_iact_stars_inject;
absorption_sum_tot += p->rt_data.debug_radiation_absorbed_tot;
/* Reset all values here in case particles won't be active next step */
p->rt_data.debug_iact_stars_inject = 0;
- iacts_with_stars_sum += p->rt_data.debug_iact_stars_inject_prep;
- iacts_with_stars_sum_tot += p->rt_data.debug_iact_stars_inject_prep_tot;
- /* Reset all values here in case particles won't be active next step */
- p->rt_data.debug_iact_stars_inject_prep = 0;
-
/* Sum up total energies for budget */
for (int g = 0; g < RT_NGROUPS; g++) {
energy_sum[g] +=
@@ -134,12 +131,9 @@ static void rt_debugging_end_of_step_hydro_mapper(void *restrict map_data,
}
}
- atomic_add(&e->rt_props->debug_radiation_absorbed_this_step, absorption_sum);
+ atomic_add(&e->rt_props->debug_radiation_absorbed_this_step,
+ absorption_sum_this_step);
atomic_add(&e->rt_props->debug_radiation_absorbed_tot, absorption_sum_tot);
- atomic_add(&e->rt_props->debug_part_injection_prep_iacts_with_stars_this_step,
- iacts_with_stars_sum);
- atomic_add(&e->rt_props->debug_part_injection_prep_iacts_with_stars_tot,
- iacts_with_stars_sum_tot);
for (int g = 0; g < RT_NGROUPS; g++) {
atomic_add_f(&(e->rt_props->debug_total_radiation_conserved_energy[g]),
@@ -160,20 +154,21 @@ rt_debugging_checks_end_of_step(struct engine *e, int verbose) {
struct space *s = e->s;
if (!(e->policy & engine_policy_rt)) return;
+#ifdef WITH_MPI
+ /* Since we aren't sending data back, none of these checks will
+ * pass a run over MPI. */
+ return;
+#endif
const ticks tic = getticks();
- /* reset values before the particle loops */
- e->rt_props->debug_radiation_emitted_this_step = 0;
- e->rt_props->debug_radiation_absorbed_this_step = 0;
- e->rt_props->debug_star_injection_prep_iacts_with_parts_this_step = 0;
- e->rt_props->debug_part_injection_prep_iacts_with_stars_this_step = 0;
- /* reset total counts as well. We track the totals since the beginning
+ /* reset values before the particle loops.
+ * reset total counts as well. We track the totals since the beginning
* of time in particles individually. */
- e->rt_props->debug_radiation_emitted_tot = 0LL;
- e->rt_props->debug_radiation_absorbed_tot = 0LL;
- e->rt_props->debug_star_injection_prep_iacts_with_parts_tot = 0LL;
- e->rt_props->debug_part_injection_prep_iacts_with_stars_tot = 0LL;
+ e->rt_props->debug_radiation_emitted_this_step = 0ULL;
+ e->rt_props->debug_radiation_absorbed_this_step = 0ULL;
+ e->rt_props->debug_radiation_emitted_tot = 0ULL;
+ e->rt_props->debug_radiation_absorbed_tot = 0ULL;
for (int g = 0; g < RT_NGROUPS; g++) {
e->rt_props->debug_total_radiation_conserved_energy[g] = 0.f;
e->rt_props->debug_total_star_emitted_energy[g] = 0.f;
@@ -191,28 +186,6 @@ rt_debugging_checks_end_of_step(struct engine *e, int verbose) {
s->sparts, s->nr_sparts, sizeof(struct spart),
threadpool_auto_chunk_size, /*extra_data=*/e);
- /* cell loop */
- /* for (int i = 0; i < s->nr_cells; ++i){ */
- /* if (s->cells_top[i].nodeID == engine_rank) { */
- /* rt_debugging_check_cell(&s->cells_top[i]); */
- /* } */
- /* } */
-
- /* message("This step: %12d %12d %12d %12d", */
- /* e->rt_props->debug_radiation_emitted_this_step, */
- /* e->rt_props->debug_radiation_absorbed_this_step, */
- /* e->rt_props->debug_star_injection_prep_iacts_with_parts_this_step,
- */
- /* e->rt_props->debug_part_injection_prep_iacts_with_stars_this_step
- */
- /* ); */
- /* message("Over lifetime: %12lld %12lld %12lld %12lld", */
- /* e->rt_props->debug_radiation_emitted_tot, */
- /* e->rt_props->debug_radiation_absorbed_tot, */
- /* e->rt_props->debug_star_injection_prep_iacts_with_parts_tot, */
- /* e->rt_props->debug_part_injection_prep_iacts_with_stars_tot */
- /* ); */
-
/* Have we accidentally invented or deleted some radiation somewhere? */
if ((e->rt_props->debug_radiation_emitted_this_step !=
e->rt_props->debug_radiation_absorbed_this_step) ||
@@ -220,39 +193,13 @@ rt_debugging_checks_end_of_step(struct engine *e, int verbose) {
e->rt_props->debug_radiation_absorbed_tot))
error(
"Emitted and absorbed radiation vary.\n"
- " This step: star emission %12d; gas absorption %12d\n"
+ " This step: star emission %12lld; gas absorption %12lld\n"
"Since start: star emission %12lld; gas absorption %12lld",
e->rt_props->debug_radiation_emitted_this_step,
e->rt_props->debug_radiation_absorbed_this_step,
e->rt_props->debug_radiation_emitted_tot,
e->rt_props->debug_radiation_absorbed_tot);
- if ((e->rt_props->debug_part_injection_prep_iacts_with_stars_this_step !=
- e->rt_props->debug_star_injection_prep_iacts_with_parts_this_step) ||
- (e->rt_props->debug_part_injection_prep_iacts_with_stars_tot !=
- e->rt_props->debug_star_injection_prep_iacts_with_parts_tot))
- error(
- "Injection prep counts parts vs stars disagree.\n"
- " This step: star iacts: %12d; gas iacts: %12d\n"
- "Since start: star iacts: %12lld; gas iacts: %12lld",
- e->rt_props->debug_star_injection_prep_iacts_with_parts_this_step,
- e->rt_props->debug_part_injection_prep_iacts_with_stars_this_step,
- e->rt_props->debug_star_injection_prep_iacts_with_parts_tot,
- e->rt_props->debug_part_injection_prep_iacts_with_stars_tot);
-
- if ((e->rt_props->debug_part_injection_prep_iacts_with_stars_this_step !=
- e->rt_props->debug_radiation_emitted_this_step) ||
- (e->rt_props->debug_part_injection_prep_iacts_with_stars_tot !=
- e->rt_props->debug_radiation_emitted_tot))
- error(
- "Injection prep iact counts vs actual iact counts disagree.\n"
- " This step: prep iacts: %12d; inject iacts: %12d\n"
- "Since start: prep iacts: %12lld; inject iacts: %12lld",
- e->rt_props->debug_part_injection_prep_iacts_with_stars_this_step,
- e->rt_props->debug_radiation_emitted_this_step,
- e->rt_props->debug_part_injection_prep_iacts_with_stars_tot,
- e->rt_props->debug_radiation_emitted_tot);
-
/* Write down energy budget for this timestep. */
if (e->step > 1) {
fprintf(e->rt_props->conserved_energy_filep, "\n");
@@ -282,42 +229,5 @@ rt_debugging_checks_end_of_step(struct engine *e, int verbose) {
clocks_getunit());
}
-/**
- * @brief This function is intended for debugging purposes only. It is called
- * during the self injection tasks, (regardless whether the particle actually
- * has neighbours to interact with) and intended to mark star or gas particles
- * to have been called during the step so further checks can be performed
- * further down the task system.
- *
- * @param p Hydro particle.
- * @param rt_props RT properties struct
- */
-__attribute__((always_inline)) INLINE static void
-rt_debugging_check_injection_part(struct part *restrict p,
- struct rt_props *props) {
-
- /* skip this for GEAR */
- /* if (props->do_all_parts_have_stars_checks) */
- /* p->rt_data.debug_injection_check += 1; */
-}
-
-/**
- * @brief This function is intended for debugging purposes only. It is called
- * during the self injection tasks, (regardless whether the particle actually
- * has neighbours to interact with) and intended to mark star or gas particles
- * to have been called during the step so further checks can be performed
- * further down the task system.
- *
- * @param s Star particle.
- * @param rt_props RT properties struct
- */
-__attribute__((always_inline)) INLINE static void
-rt_debugging_check_injection_spart(struct spart *restrict s,
- struct rt_props *props) {
-
- /* skip this for GEAR */
- /* if (props->do_all_parts_have_stars_checks) */
- /* s->rt_data.debug_injection_check += 1; */
-}
#endif /* SWIFT_RT_DEBUG_CHECKS */
#endif /* SWIFT_RT_DEBUGGING_GEAR_H */
diff --git a/src/rt/GEAR/rt_iact.h b/src/rt/GEAR/rt_iact.h
index 9d9d281f32d0242fed93dcc9912bfce18f4cc6a4..430bddd88e6d32ea652991ce60bd8d298de9059c 100644
--- a/src/rt/GEAR/rt_iact.h
+++ b/src/rt/GEAR/rt_iact.h
@@ -45,24 +45,16 @@
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_rt_injection_prep(const float r2, const float *dx,
const float hi, const float hj,
- struct spart *si, struct part *pj,
+ struct spart *si, const struct part *pj,
const struct cosmology *cosmo,
const struct rt_props *rt_props) {
- /* NOTE: `struct part *pj` should be `const struct part *pj`,
- * but I allow changes to it for debugging routines at the moment.
- * Nevertheless, you shouldn't be changing anything in a particle
- * in this function. */
-
/* If the star doesn't have any neighbours, we
* have nothing to do here. */
if (si->density.wcount == 0.f) return;
#ifdef SWIFT_RT_DEBUG_CHECKS
si->rt_data.debug_iact_hydro_inject_prep += 1;
- si->rt_data.debug_iact_hydro_inject_prep_tot += 1ULL;
- pj->rt_data.debug_iact_stars_inject_prep += 1;
- pj->rt_data.debug_iact_stars_inject_prep_tot += 1ULL;
#endif
/* Compute the weight of the neighbouring particle */
@@ -73,7 +65,11 @@ runner_iact_nonsym_rt_injection_prep(const float r2, const float *dx,
kernel_eval(xi, &wi);
const float hi_inv_dim = pow_dimension(hi_inv);
/* psi(x_star - x_gas, h_star) */
- const float psi = wi * hi_inv_dim / si->density.wcount;
+ /* Note: skip the devision by si->density.wcount here. It'll cancel out by the
+ * normalization anyway, and furthermore now that the injection prep is done
+ * during the star density loop, si->density.wcount won't be computed at this
+ * stage yet. */
+ const float psi = wi * hi_inv_dim;
/* Now add that weight to the appropriate octant */
int octant_index = 0;
@@ -111,17 +107,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_rt_inject(
* before other potential early exits */
if (si->rt_data.debug_iact_hydro_inject_prep == 0)
error(
- "Injecting energy from star that wasn't called"
- " during injection prep");
- if (pj->rt_data.debug_iact_stars_inject_prep == 0) {
+ "Injecting energy from star that wasn't called during injection prep");
- const float hig2 = hi * hi * kernel_gamma2;
- const float res = sqrtf(r2 / hig2);
- error(
- "Injecting energy into part that wasn't called"
- " during injection prep: sID %lld pID %lld r/H_s %.6f",
- si->id, pj->id, res);
- }
+ if (!si->rt_data.debug_emission_rate_set)
+ error("Injecting energy from star without setting emission rate");
si->rt_data.debug_iact_hydro_inject += 1;
si->rt_data.debug_radiation_emitted_tot += 1ULL;
@@ -129,18 +118,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_rt_inject(
pj->rt_data.debug_iact_stars_inject += 1;
pj->rt_data.debug_radiation_absorbed_tot += 1ULL;
- /* Attempt to catch race condition/dependency error */
- if (si->rt_data.debug_iact_hydro_inject_prep <
- si->rt_data.debug_iact_hydro_inject)
- error(
- "Star interacts with more particles during"
- " injection than during injection prep");
-
- if (pj->rt_data.debug_iact_stars_inject_prep <
- pj->rt_data.debug_iact_stars_inject)
- error(
- "Part interacts with more stars during"
- " injection than during injection prep");
#endif
/* Compute the weight of the neighbouring particle */
@@ -151,7 +128,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_rt_inject(
kernel_eval(xi, &wi);
const float hi_inv_dim = pow_dimension(hi_inv);
/* psi(x_star - x_gas, h_star) */
- const float psi = wi * hi_inv_dim / si->density.wcount;
+ /* Skip the division by si->density.wcount to remain consistent */
+ const float psi = wi * hi_inv_dim;
#if defined(HYDRO_DIMENSION_3D)
const int maxind = 8;
@@ -185,7 +163,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_rt_inject(
const float Vinv = 1.f / pj->geometry.volume;
/* Nurse, the patient is ready now */
- /* TODO: this is done differently for RT_HYDRO_CONTROLLED_INJECTION */
for (int g = 0; g < RT_NGROUPS; g++) {
/* Inject energy. */
const float injected_energy_density =
@@ -210,10 +187,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_rt_inject(
"Injecting abnormal energy spart %lld part %lld group %d | %.6e %.6e "
"%.6e",
si->id, pj->id, g, injected_energy, weight,
+
si->rt_data.emission_this_step[g]);
si->rt_data.debug_injected_energy[g] += injected_energy;
si->rt_data.debug_injected_energy_tot[g] += injected_energy;
}
+ si->rt_data.debug_psi_sum += psi;
#endif
}
@@ -258,6 +237,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_rt_flux_common(
if (mode == 1) {
+ if (pj->rt_data.debug_kicked != 1)
+ error("Trying to iact transport with unkicked particle %lld (count=%d)",
+ pj->id, pj->rt_data.debug_kicked);
+
if (pj->rt_data.debug_injection_done != 1)
error(
"Trying to do iact transport when "
@@ -270,9 +253,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_rt_flux_common(
"rt_finalise_gradient count is %d",
pj->rt_data.debug_gradients_done);
- if (pj->rt_data.debug_kicked != 1)
- error("Trying to iact transport with unkicked particle %lld (count=%d)",
- pj->id, pj->rt_data.debug_kicked);
pj->rt_data.debug_calls_iact_transport_interaction += 1;
}
#endif
diff --git a/src/rt/GEAR/rt_io.h b/src/rt/GEAR/rt_io.h
index 4158395eea577315753e9fd9bd0ad0889a0416a0..0d7f3d92ad298ecf92ca8102ad53ab2d88f817b2 100644
--- a/src/rt/GEAR/rt_io.h
+++ b/src/rt/GEAR/rt_io.h
@@ -175,7 +175,7 @@ INLINE static int rt_write_particles(const struct part* parts,
"Mass fractions of all constituent species");
#ifdef SWIFT_RT_DEBUG_CHECKS
- num_elements += 8;
+ num_elements += 7;
list[3] =
io_make_output_field("RTDebugInjectionDone", INT, 1, UNIT_CONV_NO_UNITS,
0, parts, rt_data.debug_injection_done,
@@ -207,11 +207,6 @@ INLINE static int rt_write_particles(const struct part* parts,
"RTDebugRadAbsorbedTot", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0, parts,
rt_data.debug_radiation_absorbed_tot,
"Radiation absorbed by this part during its lifetime");
- list[10] = io_make_output_field("RTDebugStarsInjectPrepTotCounts", ULONGLONG,
- 1, UNIT_CONV_NO_UNITS, 0, parts,
- rt_data.debug_iact_stars_inject_prep_tot,
- "Total interactions with stars during "
- "injection prep during its lifetime");
#endif
return num_elements;
@@ -232,22 +227,21 @@ INLINE static int rt_write_stars(const struct spart* sparts,
#ifdef SWIFT_RT_DEBUG_CHECKS
num_elements += 4;
- list[0] = io_make_output_field(
+ list[0] = io_make_output_field("RTDebugHydroIact", INT, 1, UNIT_CONV_NO_UNITS,
+ 0, sparts, rt_data.debug_iact_hydro_inject,
+ "number of interactions between this star "
+ "particle and any particle during injection");
+ list[1] = io_make_output_field(
"RTDebugEmissionRateSet", INT, 1, UNIT_CONV_NO_UNITS, 0, sparts,
rt_data.debug_emission_rate_set, "Stellar photon emission rates set?");
- list[1] = io_make_output_field(
+ list[2] = io_make_output_field(
"RTDebugRadEmittedTot", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0, sparts,
rt_data.debug_radiation_emitted_tot,
"Total radiation emitted during the lifetime of this star");
- list[2] = io_make_output_field("RTDebugInjectedPhotonEnergy", FLOAT,
+ list[3] = io_make_output_field("RTDebugInjectedPhotonEnergy", FLOAT,
RT_NGROUPS, UNIT_CONV_ENERGY, 0, sparts,
rt_data.debug_injected_energy_tot,
"Total radiation actually injected into gas");
- list[3] = io_make_output_field("RTDebugHydroInjectPrepCountsTot", ULONGLONG,
- 1, UNIT_CONV_NO_UNITS, 0, sparts,
- rt_data.debug_iact_hydro_inject_prep_tot,
- "Total interactions with particles during "
- "injection prep during its lifetime");
#endif
return num_elements;
@@ -273,12 +267,7 @@ INLINE static void rt_write_flavour(hid_t h_grp, hid_t h_grp_columns,
/* Write scheme name */
/* ----------------- */
- if (rtp->hydro_controlled_injection) {
- io_write_attribute_s(h_grp, "RT Scheme",
- RT_IMPLEMENTATION ", hydro controlled injection");
- } else {
- io_write_attribute_s(h_grp, "RT Scheme", RT_IMPLEMENTATION);
- }
+ io_write_attribute_s(h_grp, "RT Scheme", RT_IMPLEMENTATION);
/* Write photon group counts */
/* ------------------------- */
diff --git a/src/rt/GEAR/rt_properties.h b/src/rt/GEAR/rt_properties.h
index 8d86d157a18868107a2496c5d0e7fedca3429f0b..66cf1b6f090ee593eb525d74ac4fe9a83234b9cc 100644
--- a/src/rt/GEAR/rt_properties.h
+++ b/src/rt/GEAR/rt_properties.h
@@ -48,15 +48,6 @@ static void rt_interaction_rates_init(struct rt_props* restrict rt_props,
*/
struct rt_props {
- /* Are we running with hydro or star controlled injection?
- * This is added to avoid #ifdef macros as far as possible */
- int hydro_controlled_injection;
-
- /* Do we need to run a conversion after the zeroth
- * step, but before the first step? */
- int convert_stars_after_zeroth_step;
- int convert_parts_after_zeroth_step;
-
/* Are we using constant stellar emission rates? */
int use_const_emission_rates;
@@ -119,15 +110,10 @@ struct rt_props {
chemistry_data_storage* grackle_chemistry_rates;
#ifdef SWIFT_RT_DEBUG_CHECKS
- /* Do extended tests where we assume that all parts
- * have spart neighbours? */
- /* skip this for GEAR */
- /* int debug_do_all_parts_have_stars_checks; */
-
/* radiation emitted by stars this step. This is not really a property,
* but a placeholder to sum up a global variable. It's being reset
* every timestep. */
- int debug_radiation_emitted_this_step;
+ unsigned long long debug_radiation_emitted_this_step;
/* total radiation emitted by stars. This is not really a property,
* but a placeholder to sum up a global variable */
@@ -135,31 +121,14 @@ struct rt_props {
/* radiation absorbed by gas this step. This is not really a property,
* but a placeholder to sum up a global variable */
- int debug_radiation_absorbed_this_step;
+ unsigned long long debug_radiation_absorbed_this_step;
/* total radiation absorbed by gas. This is not really a property,
* but a placeholder to sum up a global variable */
unsigned long long debug_radiation_absorbed_tot;
- /* Interactions of a star with gas during injection prep this step. This is
- * not really a property, but a placeholder to sum up a global variable */
- int debug_star_injection_prep_iacts_with_parts_this_step;
-
- /* Interactions of a star with gas during injection prep. This is not
- * really a property, but a placeholder to sum up a global variable */
- unsigned long long debug_star_injection_prep_iacts_with_parts_tot;
-
- /* Interactions of a star with gas during injection prep this step. This is
- * not really a property, but a placeholder to sum up a global variable */
- int debug_part_injection_prep_iacts_with_stars_this_step;
-
- /* Interactions of a star with gas during injection prep. This is not
- * really a property, but a placeholder to sum up a global variable */
- unsigned long long debug_part_injection_prep_iacts_with_stars_tot;
-
/* Total radiation energy in the gas. It's being reset every step. */
float debug_total_radiation_conserved_energy[RT_NGROUPS];
- float debug_total_radiation_energy_density[RT_NGROUPS];
float debug_total_star_emitted_energy[RT_NGROUPS];
/* Files to write energy budget to after every step */
@@ -168,6 +137,45 @@ struct rt_props {
#endif
};
+/**
+ * @brief open up files to write some debugging check outputs.
+ * This function is temporary for development, and shouldn't stay
+ * long.
+ *
+ * @param rtp #rt_props struct
+ * @param mode open files with this mode. "w" for new file, "a" for append.
+ **/
+static void rt_props_open_debugging_files(struct rt_props* rtp,
+ const char* mode) {
+#ifdef SWIFT_RT_DEBUG_CHECKS
+#ifdef WITH_MPI
+ return;
+#endif
+
+ rtp->conserved_energy_filep = fopen("RT_conserved_energy_budget.txt", mode);
+ if (rtp->conserved_energy_filep == NULL)
+ error("Couldn't open RT conserved energy budget file to write in");
+ rtp->star_emitted_energy_filep = fopen("RT_star_injected_energy.txt", mode);
+ if (rtp->star_emitted_energy_filep == NULL)
+ error("Couldn't open RT star energy budget file to write in");
+
+ if (strcmp(mode, "w") == 0 && rtp->use_const_emission_rates) {
+ /* If we're starting a new file, dump the header first */
+ FILE* files[2] = {rtp->conserved_energy_filep,
+ rtp->star_emitted_energy_filep};
+ for (int f = 0; f < 2; f++) {
+ fprintf(files[f], "# Emission rates: ");
+ const double solar_luminosity = 3.826e33; /* erg/s */
+ for (int g = 0; g < RT_NGROUPS; g++) {
+ fprintf(files[f], "%12.6e ",
+ rtp->stellar_const_emission_rates[g] * solar_luminosity);
+ }
+ fprintf(files[f], "\n");
+ }
+ }
+#endif
+};
+
/**
* @brief initialize grackle during rt_props_init
*
@@ -292,21 +300,8 @@ __attribute__((always_inline)) INLINE static void rt_props_init(
const struct unit_system* us, struct swift_params* params,
struct cosmology* cosmo) {
-#ifdef RT_HYDRO_CONTROLLED_INJECTION
- rtp->hydro_controlled_injection = 1;
-#else
- rtp->hydro_controlled_injection = 0;
-#endif
-
/* Make sure we reset debugging counters correctly after
* zeroth step. */
-#ifdef SWIFT_RT_DEBUG_CHECKS
- rtp->convert_parts_after_zeroth_step = 1;
- rtp->convert_stars_after_zeroth_step = 1;
-#else
- rtp->convert_parts_after_zeroth_step = 0;
- rtp->convert_stars_after_zeroth_step = rtp->hydro_controlled_injection;
-#endif
/* Read in photon frequency group properties */
/* ----------------------------------------- */
@@ -468,32 +463,16 @@ __attribute__((always_inline)) INLINE static void rt_props_init(
#ifdef SWIFT_RT_DEBUG_CHECKS
rtp->debug_radiation_emitted_tot = 0ULL;
+ rtp->debug_radiation_emitted_this_step = 0ULL;
+
rtp->debug_radiation_absorbed_tot = 0ULL;
- rtp->debug_star_injection_prep_iacts_with_parts_tot = 0LL;
- rtp->debug_part_injection_prep_iacts_with_stars_tot = 0LL;
+ rtp->debug_radiation_absorbed_this_step = 0ULL;
+
for (int g = 0; g < RT_NGROUPS; g++)
rtp->debug_total_star_emitted_energy[g] = 0.f;
- /* Open up files for energy budgets */
- rtp->conserved_energy_filep = fopen("RT_conserved_energy_budget.txt", "w");
- if (rtp->conserved_energy_filep == NULL)
- error("Couldn't open RT conserved energy budget file to write in");
- rtp->star_emitted_energy_filep = fopen("RT_star_injected_energy.txt", "w");
- if (rtp->star_emitted_energy_filep == NULL)
- error("Couldn't open RT star energy budget file to write in");
+ rt_props_open_debugging_files(rtp, "w");
- if (rtp->use_const_emission_rates) {
- FILE* files[2] = {rtp->conserved_energy_filep,
- rtp->star_emitted_energy_filep};
- for (int f = 0; f < 2; f++) {
- fprintf(files[f], "# Emission rates: ");
- const double solar_luminosity = 3.826e33; /* erg/s */
- for (int g = 0; g < RT_NGROUPS; g++)
- fprintf(files[f], "%12.6e ",
- rtp->stellar_const_emission_rates[g] * solar_luminosity);
- fprintf(files[f], "\n");
- }
- }
#endif
/* Grackle setup */
@@ -528,6 +507,10 @@ __attribute__((always_inline)) INLINE static void rt_struct_dump(
restart_write_blocks((void*)props, sizeof(struct rt_props), 1, stream,
"RT props", "RT properties struct");
+ /* The RT parameters, in particular the reduced speed of light, are
+ * not defined at compile time. So we need to read them in again. */
+ restart_write_blocks(&rt_params, sizeof(struct rt_parameters), 1, stream,
+ "RT global parameters", "RT global parameters struct");
}
/**
@@ -542,6 +525,14 @@ __attribute__((always_inline)) INLINE static void rt_struct_restore(
restart_read_blocks((void*)props, sizeof(struct rt_props), 1, stream, NULL,
"RT properties struct");
+ /* The RT parameters, in particular the reduced speed of light, are
+ * not defined at compile time. So we need to write them down. */
+ restart_read_blocks(&rt_params, sizeof(struct rt_parameters), 1, stream, NULL,
+ "RT global parameters struct");
+#ifdef SWIFT_RT_DEBUG_CHECKS
+ /* Reset the file pointers for temporary stats */
+ rt_props_open_debugging_files(props, "a");
+#endif
}
#endif /* SWIFT_RT_PROPERTIES_GEAR_H */
diff --git a/src/rt/GEAR/rt_stellar_emission_rate.h b/src/rt/GEAR/rt_stellar_emission_rate.h
index a6febb1df0fead661d08bb15a6e6961baee484f8..f9980925fe7dc7f3f83c4a8fc9f767b7e8f353a8 100644
--- a/src/rt/GEAR/rt_stellar_emission_rate.h
+++ b/src/rt/GEAR/rt_stellar_emission_rate.h
@@ -55,11 +55,9 @@ __attribute__((always_inline)) INLINE static void rt_set_stellar_emission_rate(
rt_props->stellar_const_emission_rates[g] * solar_luminosity;
sp->rt_data.emission_this_step[g] = emission_rate_internal_units * dt;
}
+ } else {
+ error("Unknown stellar emission rate model");
}
-
-#ifdef SWIFT_RT_DEBUG_CHECKS
- sp->rt_data.debug_emission_rate_set += 1;
-#endif
}
#endif /* SWIFT_RT_STELLAR_EMISSION_RATE_GEAR_H */
diff --git a/src/rt/GEAR/rt_struct.h b/src/rt/GEAR/rt_struct.h
index bd3f6220459a2293c740568f81b0c86be918523f..b02d9f5f6054105fed150952acf5ff91e77901d2 100644
--- a/src/rt/GEAR/rt_struct.h
+++ b/src/rt/GEAR/rt_struct.h
@@ -85,26 +85,13 @@ struct rt_part_data {
/*! how much radiation this part received from stars during total lifetime */
unsigned long long debug_radiation_absorbed_tot;
- /*! how many interactions this part had with stars in injection prep over
- * total lifetime */
- unsigned long long debug_iact_stars_inject_prep_tot;
-
/* data to store during one time step */
- /*! how many stars this part interacted with during preparation*/
- /* Note: It's useless to write this in outputs, as it gets reset
- * at the end of every step. */
- int debug_iact_stars_inject_prep;
-
/*! how many stars this part interacted with during injection*/
/* Note: It's useless to write this in outputs, as it gets reset
* at the end of every step. */
int debug_iact_stars_inject;
- /* skip this for GEAR */
- /* called in a self/rt_injection task? */
- /* int debug_injection_check; */
-
/*! calls from gradient interaction loop in actual function */
int debug_calls_iact_gradient_interaction;
@@ -136,8 +123,6 @@ struct rt_spart_data {
/* Stellar energy emission that will be injected in to gas.
* Total energy, not density, not rate! */
- /* TODO: keep this also for RT_HYDRO_CONTROLLED_INJECTION and
- * store results with each hydro-star interaction in here */
float emission_this_step[RT_NGROUPS];
/*! Neighbour weigths in each octant surrounding the star */
@@ -149,36 +134,28 @@ struct rt_spart_data {
/*! how much radiation this star emitted during total lifetime */
unsigned long long debug_radiation_emitted_tot;
- /*! how many interactions this star had with parts during
- * injection prep over total lifetime */
- unsigned long long debug_iact_hydro_inject_prep_tot;
-
/* data to store during one time step */
- /*! how many hydro particles this particle interacted with */
- /* Note: It's useless to write this in outputs, as it gets reset
- * at the end of every step. */
+ /*! how many hydro particles this particle interacted with
+ * during injection */
int debug_iact_hydro_inject;
/*! how many hydro particles this particle interacted with
* during injection prep*/
- /* Note: It's useless to write this in outputs, as it gets reset
- * at the end of every step. */
int debug_iact_hydro_inject_prep;
/*! stellar photon emisison rate computed? */
int debug_emission_rate_set;
- /* skip this for GEAR */
- /* !called in a self/rt_injection task? */
- /* int debug_injection_check; */
-
/*! how much energy this star particle actually has injected into the gas */
float debug_injected_energy[RT_NGROUPS];
/*! how much energy this star particle actually has injected into the gas over
* the entire run*/
float debug_injected_energy_tot[RT_NGROUPS];
+
+ /*! sum up total weights used during injection to compare consistency */
+ float debug_psi_sum;
#endif
};
diff --git a/src/rt/GEAR/rt_thermochemistry.h b/src/rt/GEAR/rt_thermochemistry.h
index fcb65f7191ff00815ec0b888a5ad53472d187558..55d6e7d9685579a8529aee18e2b80adf806048d5 100644
--- a/src/rt/GEAR/rt_thermochemistry.h
+++ b/src/rt/GEAR/rt_thermochemistry.h
@@ -91,20 +91,6 @@ static void rt_do_thermochemistry(struct part* restrict p,
/* Note: Can't pass rt_props as const struct because of grackle
* accessinging its properties there */
-#ifdef SWIFT_RT_DEBUG_CHECKS
- if (p->rt_data.debug_kicked != 1)
- error("Trying to do thermochemistry on unkicked particle %lld (count=%d)",
- p->id, p->rt_data.debug_kicked);
- if (!p->rt_data.debug_injection_done)
- error("Trying to do thermochemistry when injection step hasn't been done");
- if (!p->rt_data.debug_gradients_done)
- error("Trying to do thermochemistry when gradient step hasn't been done");
- if (!p->rt_data.debug_transport_done)
- error("Trying to do thermochemistry when transport step hasn't been done");
-
- p->rt_data.debug_thermochem_done += 1;
-#endif
-
/* Nothing to do here? */
if (rt_props->skip_thermochemistry) return;
if (dt == 0.) return;
diff --git a/src/rt/SPHM1RT/rt.h b/src/rt/SPHM1RT/rt.h
index b4e0d0922b9e6aa17abee3ad5501a47f5135a6d2..fa22669745259f9afa5fd40c7bea31fef6acad96 100644
--- a/src/rt/SPHM1RT/rt.h
+++ b/src/rt/SPHM1RT/rt.h
@@ -285,7 +285,7 @@ __attribute__((always_inline)) INLINE static void rt_first_init_part(
/**
* @brief Initialises particle quantities that can't be set
* otherwise before the zeroth step is finished. E.g. because
- * they require the particle density to be known.
+ * they require the particle density and time step to be known.
*
* @param p particle to work on
* @param rt_props RT properties struct
@@ -321,7 +321,7 @@ __attribute__((always_inline)) INLINE static void rt_first_init_spart(
/**
* @brief Initialises particle quantities that can't be set
* otherwise before the zeroth step is finished. E.g. because
- * they require the star density to be known.
+ * they require the star density and time step to be known.
* @param sp star particle to work on
* @param time current system time
* @param star_age age of the star *at the end of the step*
@@ -685,8 +685,9 @@ __attribute__((always_inline)) INLINE static void rt_prepare_force(
* @brief Clean the allocated memory inside the RT properties struct.
*
* @param props the #rt_props.
+ * @param restart did we restart?
*/
__attribute__((always_inline)) INLINE static void rt_clean(
- struct rt_props* props) {}
+ struct rt_props* props, int restart) {}
#endif /* SWIFT_RT_SPHM1RT_H */
diff --git a/src/rt/SPHM1RT/rt_properties.h b/src/rt/SPHM1RT/rt_properties.h
index 04909a204f331338cb1074bcf42cb9253a6839af..4f297a6ea12a4e854feddb47576442303cf5c1c3 100644
--- a/src/rt/SPHM1RT/rt_properties.h
+++ b/src/rt/SPHM1RT/rt_properties.h
@@ -46,15 +46,6 @@ struct rt_props {
* Includes 0 as leftmost edge, doesn't include infinity as
* rightmost bin edge*/
float photon_groups[RT_NGROUPS];
-
- /* Are we running with hydro or star controlled injection?
- * This is added to avoid #ifdef macros as far as possible */
- int hydro_controlled_injection;
-
- /* Do we need to run a conversion after the zeroth
- * step, but before the first step? */
- int convert_stars_after_zeroth_step;
- int convert_parts_after_zeroth_step;
};
/**
@@ -85,15 +76,6 @@ __attribute__((always_inline)) INLINE static void rt_props_init(
const struct unit_system* us, struct swift_params* params,
struct cosmology* cosmo) {
-#ifdef RT_HYDRO_CONTROLLED_INJECTION
- rtp->hydro_controlled_injection = 1;
-#else
- rtp->hydro_controlled_injection = 0;
-#endif
-
- rtp->convert_parts_after_zeroth_step = 0;
- rtp->convert_stars_after_zeroth_step = 0;
-
if (RT_NGROUPS <= 0) {
error(
"You need to run SPHM1RT with at least 1 photon group, "
diff --git a/src/rt/debug/rt.h b/src/rt/debug/rt.h
index b5115ab5fb2ec1de0dc0f165aa3c212342a01f9a..1a1ebb7d12b374960364296a696726184f1b8d65 100644
--- a/src/rt/debug/rt.h
+++ b/src/rt/debug/rt.h
@@ -20,8 +20,6 @@
#define SWIFT_RT_DEBUG_H
#include "rt_debugging.h"
-#include "rt_stellar_emission_rate.h"
-#include "rt_thermochemistry.h"
/**
* @file src/rt/debug/rt.h
@@ -29,7 +27,7 @@
*/
/**
- * @brief Compute the photon emission rates for this stellar particle
+ * @brief Compute the photon emission rates for this stellar particle.
* This function is called every time the spart is being reset
* (during start-up and during stars ghost if spart is active)
* and assumes that the photon emission rate is an intrinsic
@@ -50,20 +48,11 @@ rt_compute_stellar_emission_rate(struct spart* restrict sp, double time,
const struct phys_const* phys_const,
const struct unit_system* internal_units) {
- /* Skip initial fake time-step */
- if (dt == 0.0l) return;
+ sp->rt_data.debug_emission_rate_set += 1;
- if (time == 0.l) {
- /* if function is called before the first actual step, time is still
- * at zero unless specified otherwise in parameter file.*/
- star_age = dt;
- }
-
- /* now get the emission rates */
- double star_age_begin_of_step = star_age - dt;
- star_age_begin_of_step = max(0.l, star_age_begin_of_step);
- rt_set_stellar_emission_rate(sp, star_age_begin_of_step, star_age, rt_props,
- phys_const, internal_units);
+ /* rt_set_stellar_emission_rate(sp, star_age_begin_of_step, star_age,
+ * rt_props, */
+ /* phys_const, internal_units); */
}
/**
@@ -90,9 +79,7 @@ __attribute__((always_inline)) INLINE static void rt_reset_part(
/* reset this here as well as in the rt_debugging_checks_end_of_step()
* routine to test task dependencies are done right */
p->rt_data.debug_iact_stars_inject = 0;
- p->rt_data.debug_iact_stars_inject_prep = 0;
- p->rt_data.debug_injection_check = 0;
p->rt_data.debug_calls_iact_gradient_interaction = 0;
p->rt_data.debug_calls_iact_transport_interaction = 0;
@@ -114,13 +101,12 @@ __attribute__((always_inline)) INLINE static void rt_first_init_part(
rt_init_part(p);
rt_reset_part(p);
p->rt_data.debug_radiation_absorbed_tot = 0ULL;
- p->rt_data.debug_iact_stars_inject_prep_tot = 0ULL;
}
/**
* @brief Initialises particle quantities that can't be set
* otherwise before the zeroth step is finished. E.g. because
- * they require the particle density to be known.
+ * they require the particle density and time step to be known.
*
* @param p particle to work on
* @param rt_props RT properties struct
@@ -132,8 +118,16 @@ rt_init_part_after_zeroth_step(struct part* restrict p,
/* If we're running with debugging checks on, reset debugging
* counters and flags in particular after the zeroth step so
* that the checks work as intended. */
+ rt_init_part(p);
rt_reset_part(p);
+ /* Since the inject_prep has been moved to the density loop, the
+ * initialization at startup is messing with the total counters for stars
+ * because the density is called, but not the force-and-kick tasks. So reset
+ * the total counters here as well so that they will match the star counters.
+ */
+ p->rt_data.debug_radiation_absorbed_tot = 0ULL;
}
+
/**
* @brief Initialisation of the RT density loop related star particle data.
* Note: during initalisation (space_init), rt_reset_spart and rt_init_spart
@@ -142,7 +136,14 @@ rt_init_part_after_zeroth_step(struct part* restrict p,
* @param sp star particle to work on
*/
__attribute__((always_inline)) INLINE static void rt_init_spart(
- struct spart* restrict sp) {}
+ struct spart* restrict sp) {
+
+ /* reset this here as well as in the rt_debugging_checks_end_of_step()
+ * routine to test task dependencies are done right */
+ sp->rt_data.debug_iact_hydro_inject_prep = 0;
+ sp->rt_data.debug_iact_hydro_inject = 0;
+ sp->rt_data.debug_emission_rate_set = 0;
+}
/**
* @brief Reset of the RT star particle data not related to the density.
@@ -153,18 +154,7 @@ __attribute__((always_inline)) INLINE static void rt_init_spart(
* @param sp star particle to work on
*/
__attribute__((always_inline)) INLINE static void rt_reset_spart(
- struct spart* restrict sp) {
-
- /* reset everything */
-
- /* reset this here as well as in the rt_debugging_checks_end_of_step()
- * routine to test task dependencies are done right */
- sp->rt_data.debug_iact_hydro_inject = 0;
- sp->rt_data.debug_iact_hydro_inject_prep = 0;
-
- sp->rt_data.debug_emission_rate_set = 0;
- sp->rt_data.debug_injection_check = 0;
-}
+ struct spart* restrict sp) {}
/**
* @brief First initialisation of the RT star particle data.
@@ -182,7 +172,7 @@ __attribute__((always_inline)) INLINE static void rt_first_init_spart(
/**
* @brief Initialises particle quantities that can't be set
* otherwise before the zeroth step is finished. E.g. because
- * they require the star density to be known.
+ * they require the star density and time step to be known.
* @param sp star particle to work on
* @param time current system time
* @param star_age age of the star *at the end of the step*
@@ -201,7 +191,13 @@ rt_init_star_after_zeroth_step(struct spart* restrict sp, double time,
/* If we're running with debugging checks on, reset debugging
* counters and flags in particular after the zeroth step so
* that the checks work as intended. */
+ rt_init_spart(sp);
rt_reset_spart(sp);
+ /* Since the inject_prep has been moved to the density loop, the
+ * initialization at startup is messing with the total counters because
+ * the density is called, but not the force-and-kick tasks. So reset
+ * the total counters here as well. */
+ sp->rt_data.debug_radiation_emitted_tot = 0ULL;
}
/**
@@ -307,12 +303,8 @@ __attribute__((always_inline)) INLINE static double rt_part_dt(
__attribute__((always_inline)) INLINE static void rt_finalise_injection(
struct part* restrict p, struct rt_props* props) {
- if (props->debug_do_all_parts_have_stars_checks &&
- p->rt_data.debug_injection_check != 1)
- error("called ghost1 when injection check count is %d; ID=%lld",
- p->rt_data.debug_injection_check, p->id);
if (p->rt_data.debug_kicked != 1)
- error("called ghost1 when particle %lld is unkicked (count=%d)", p->id,
+ error("called rt_ghost1 when particle %lld is unkicked (count=%d)", p->id,
p->rt_data.debug_kicked);
p->rt_data.debug_injection_done += 1;
}
@@ -399,7 +391,19 @@ __attribute__((always_inline)) INLINE static void rt_tchem(
const struct phys_const* restrict phys_const,
const struct unit_system* restrict us, const double dt) {
- rt_do_thermochemistry(p);
+ if (p->rt_data.debug_kicked != 1)
+ error("Trying to do thermochemistry on unkicked particle %lld (count=%d)",
+ p->id, p->rt_data.debug_kicked);
+ if (!p->rt_data.debug_injection_done)
+ error("Trying to do thermochemistry when injection step hasn't been done");
+ if (!p->rt_data.debug_gradients_done)
+ error("Trying to do thermochemistry when gradient step hasn't been done");
+ if (!p->rt_data.debug_transport_done)
+ error("Trying to do thermochemistry when transport step hasn't been done");
+
+ p->rt_data.debug_thermochem_done += 1;
+
+ /* rt_do_thermochemistry(p); */
}
/**
@@ -429,7 +433,7 @@ __attribute__((always_inline)) INLINE static void rt_kick_extra(
/**
* @brief Prepare a particle for the !HYDRO! force calculation.
* E.g. for the meshless schemes, we need to take into account the
- * mass fluxes of the ionizing species between particles.
+ * mass fluxes of the constituent species between particles.
* NOTE: don't call this during rt_init_part or rt_reset_part,
* follow the hydro_prepare_force logic.
*
@@ -442,8 +446,9 @@ __attribute__((always_inline)) INLINE static void rt_prepare_force(
* @brief Clean the allocated memory inside the RT properties struct.
*
* @param props the #rt_props.
+ * @param restart did we restart?
*/
__attribute__((always_inline)) INLINE static void rt_clean(
- struct rt_props* props) {}
+ struct rt_props* props, int restart) {}
#endif /* SWIFT_RT_DEBUG_H */
diff --git a/src/rt/debug/rt_debugging.h b/src/rt/debug/rt_debugging.h
index ae52fda0e08ec1bae62bd9520eadc03ceeae3fe7..83e670d3f22f8a018e76c654751f95a701c55a85 100644
--- a/src/rt/debug/rt_debugging.h
+++ b/src/rt/debug/rt_debugging.h
@@ -37,31 +37,22 @@ static void rt_debugging_end_of_step_stars_mapper(void *restrict map_data,
struct spart *restrict sparts = (struct spart *)map_data;
const struct engine *restrict e = (struct engine *)extra_data;
- int emission_sum = 0;
- int iacts_with_parts_sum = 0;
+ unsigned long long emission_sum_this_step = 0ULL;
unsigned long long emission_sum_tot = 0ULL;
- unsigned long long iacts_with_parts_sum_tot = 0ULL;
for (int k = 0; k < scount; k++) {
struct spart *restrict sp = &sparts[k];
- emission_sum += sp->rt_data.debug_iact_hydro_inject;
+ emission_sum_this_step += sp->rt_data.debug_iact_hydro_inject;
emission_sum_tot += sp->rt_data.debug_radiation_emitted_tot;
/* Reset all values here in case stars won't be active next step */
sp->rt_data.debug_iact_hydro_inject = 0;
-
- iacts_with_parts_sum += sp->rt_data.debug_iact_hydro_inject_prep;
- iacts_with_parts_sum_tot += sp->rt_data.debug_iact_hydro_inject_prep_tot;
- /* Reset all values here in case stars won't be active next step */
sp->rt_data.debug_iact_hydro_inject_prep = 0;
}
- atomic_add(&e->rt_props->debug_radiation_emitted_this_step, emission_sum);
+ atomic_add(&e->rt_props->debug_radiation_emitted_this_step,
+ emission_sum_this_step);
atomic_add(&e->rt_props->debug_radiation_emitted_tot, emission_sum_tot);
- atomic_add(&e->rt_props->debug_star_injection_prep_iacts_with_parts_this_step,
- iacts_with_parts_sum);
- atomic_add(&e->rt_props->debug_star_injection_prep_iacts_with_parts_tot,
- iacts_with_parts_sum_tot);
}
/**
@@ -74,31 +65,21 @@ static void rt_debugging_end_of_step_hydro_mapper(void *restrict map_data,
struct part *restrict parts = (struct part *)map_data;
const struct engine *restrict e = (struct engine *)extra_data;
- int absorption_sum = 0;
- int iacts_with_stars_sum = 0;
+ unsigned long long absorption_sum_this_step = 0ULL;
unsigned long long absorption_sum_tot = 0ULL;
- unsigned long long iacts_with_stars_sum_tot = 0ULL;
for (int k = 0; k < count; k++) {
struct part *restrict p = &parts[k];
- absorption_sum += p->rt_data.debug_iact_stars_inject;
+ absorption_sum_this_step += p->rt_data.debug_iact_stars_inject;
absorption_sum_tot += p->rt_data.debug_radiation_absorbed_tot;
/* Reset all values here in case particles won't be active next step */
p->rt_data.debug_iact_stars_inject = 0;
-
- iacts_with_stars_sum += p->rt_data.debug_iact_stars_inject_prep;
- iacts_with_stars_sum_tot += p->rt_data.debug_iact_stars_inject_prep_tot;
- /* Reset all values here in case particles won't be active next step */
- p->rt_data.debug_iact_stars_inject_prep = 0;
}
- atomic_add(&e->rt_props->debug_radiation_absorbed_this_step, absorption_sum);
+ atomic_add(&e->rt_props->debug_radiation_absorbed_this_step,
+ absorption_sum_this_step);
atomic_add(&e->rt_props->debug_radiation_absorbed_tot, absorption_sum_tot);
- atomic_add(&e->rt_props->debug_part_injection_prep_iacts_with_stars_this_step,
- iacts_with_stars_sum);
- atomic_add(&e->rt_props->debug_part_injection_prep_iacts_with_stars_tot,
- iacts_with_stars_sum_tot);
}
/**
@@ -114,20 +95,21 @@ rt_debugging_checks_end_of_step(struct engine *e, int verbose) {
struct space *s = e->s;
if (!(e->policy & engine_policy_rt)) return;
+#ifdef WITH_MPI
+ /* Since we aren't sending data back, none of these checks will
+ * pass a run over MPI. */
+ return;
+#endif
const ticks tic = getticks();
- /* reset values before the particle loops */
- e->rt_props->debug_radiation_emitted_this_step = 0;
- e->rt_props->debug_radiation_absorbed_this_step = 0;
- e->rt_props->debug_star_injection_prep_iacts_with_parts_this_step = 0;
- e->rt_props->debug_part_injection_prep_iacts_with_stars_this_step = 0;
- /* reset total counts as well. We track the totals since the beginning
+ /* reset values before the particle loops.
+ * reset total counts as well. We track the totals since the beginning
* of time in particles individually. */
- e->rt_props->debug_radiation_emitted_tot = 0LL;
- e->rt_props->debug_radiation_absorbed_tot = 0LL;
- e->rt_props->debug_star_injection_prep_iacts_with_parts_tot = 0LL;
- e->rt_props->debug_part_injection_prep_iacts_with_stars_tot = 0LL;
+ e->rt_props->debug_radiation_emitted_this_step = 0ULL;
+ e->rt_props->debug_radiation_absorbed_this_step = 0ULL;
+ e->rt_props->debug_radiation_emitted_tot = 0ULL;
+ e->rt_props->debug_radiation_absorbed_tot = 0ULL;
/* hydro particle loop */
if (s->nr_parts > 0)
@@ -141,21 +123,6 @@ rt_debugging_checks_end_of_step(struct engine *e, int verbose) {
s->sparts, s->nr_sparts, sizeof(struct spart),
threadpool_auto_chunk_size, /*extra_data=*/e);
- /* message("This step: %12d %12d %12d %12d", */
- /* e->rt_props->debug_radiation_emitted_this_step, */
- /* e->rt_props->debug_radiation_absorbed_this_step, */
- /* e->rt_props->debug_star_injection_prep_iacts_with_parts_this_step,
- */
- /* e->rt_props->debug_part_injection_prep_iacts_with_stars_this_step
- */
- /* ); */
- /* message("Over lifetime: %12lld %12lld %12lld %12lld", */
- /* e->rt_props->debug_radiation_emitted_tot, */
- /* e->rt_props->debug_radiation_absorbed_tot, */
- /* e->rt_props->debug_star_injection_prep_iacts_with_parts_tot, */
- /* e->rt_props->debug_part_injection_prep_iacts_with_stars_tot */
- /* ); */
-
/* Have we accidentally invented or deleted some radiation somewhere? */
if ((e->rt_props->debug_radiation_emitted_this_step !=
e->rt_props->debug_radiation_absorbed_this_step) ||
@@ -163,77 +130,16 @@ rt_debugging_checks_end_of_step(struct engine *e, int verbose) {
e->rt_props->debug_radiation_absorbed_tot))
error(
"Emitted and absorbed radiation vary.\n"
- " This step: star emission %12d; gas absorption %12d\n"
+ " This step: star emission %12lld; gas absorption %12lld\n"
"Since start: star emission %12lld; gas absorption %12lld",
e->rt_props->debug_radiation_emitted_this_step,
e->rt_props->debug_radiation_absorbed_this_step,
e->rt_props->debug_radiation_emitted_tot,
e->rt_props->debug_radiation_absorbed_tot);
- if ((e->rt_props->debug_part_injection_prep_iacts_with_stars_this_step !=
- e->rt_props->debug_star_injection_prep_iacts_with_parts_this_step) ||
- (e->rt_props->debug_part_injection_prep_iacts_with_stars_tot !=
- e->rt_props->debug_star_injection_prep_iacts_with_parts_tot))
- error(
- "Injection prep counts parts vs stars disagree.\n"
- " This step: star iacts: %12d; gas iacts: %12d\n"
- "Since start: star iacts: %12lld; gas iacts: %12lld",
- e->rt_props->debug_star_injection_prep_iacts_with_parts_this_step,
- e->rt_props->debug_part_injection_prep_iacts_with_stars_this_step,
- e->rt_props->debug_star_injection_prep_iacts_with_parts_tot,
- e->rt_props->debug_part_injection_prep_iacts_with_stars_tot);
-
- if ((e->rt_props->debug_part_injection_prep_iacts_with_stars_this_step !=
- e->rt_props->debug_radiation_emitted_this_step) ||
- (e->rt_props->debug_part_injection_prep_iacts_with_stars_tot !=
- e->rt_props->debug_radiation_emitted_tot))
- error(
- "Injection prep iact counts vs actual iact counts disagree.\n"
- " This step: prep iacts: %12d; inject iacts: %12d\n"
- "Since start: prep iacts: %12lld; inject iacts: %12lld",
- e->rt_props->debug_part_injection_prep_iacts_with_stars_this_step,
- e->rt_props->debug_radiation_emitted_this_step,
- e->rt_props->debug_part_injection_prep_iacts_with_stars_tot,
- e->rt_props->debug_radiation_emitted_tot);
-
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
clocks_getunit());
}
-/**
- * @brief This function is intended for debugging purposes only. It is called
- * during the self injection tasks, (regardless whether the particle actually
- * has neighbours to interact with) and intended to mark star or gas particles
- * to have been called during the step so further checks can be performed
- * further down the task system.
- *
- * @param p Hydro particle.
- * @param rt_props RT properties struct
- */
-__attribute__((always_inline)) INLINE static void
-rt_debugging_check_injection_part(struct part *restrict p,
- struct rt_props *props) {
-
- if (props->debug_do_all_parts_have_stars_checks)
- p->rt_data.debug_injection_check += 1;
-}
-
-/**
- * @brief This function is intended for debugging purposes only. It is called
- * during the self injection tasks, (regardless whether the particle actually
- * has neighbours to interact with) and intended to mark star or gas particles
- * to have been called during the step so further checks can be performed
- * further down the task system.
- *
- * @param s Star particle.
- * @param rt_props RT properties struct
- */
-__attribute__((always_inline)) INLINE static void
-rt_debugging_check_injection_spart(struct spart *restrict s,
- struct rt_props *props) {
-
- if (props->debug_do_all_parts_have_stars_checks)
- s->rt_data.debug_injection_check += 1;
-}
#endif /* SWIFT_RT_DEBUGGING_DEBUG_H */
diff --git a/src/rt/debug/rt_iact.h b/src/rt/debug/rt_iact.h
index 4e301d2e5600b5c8edb97cf99c8c06f5164ee812..375e9804b4c643c42f0b2987e1eb2de32a908d43 100644
--- a/src/rt/debug/rt_iact.h
+++ b/src/rt/debug/rt_iact.h
@@ -44,14 +44,11 @@
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_rt_injection_prep(const float r2, const float *dx,
const float hi, const float hj,
- struct spart *si, struct part *pj,
+ struct spart *si, const struct part *pj,
const struct cosmology *cosmo,
const struct rt_props *rt_props) {
si->rt_data.debug_iact_hydro_inject_prep += 1;
- si->rt_data.debug_iact_hydro_inject_prep_tot += 1ULL;
- pj->rt_data.debug_iact_stars_inject_prep += 1;
- pj->rt_data.debug_iact_stars_inject_prep_tot += 1ULL;
}
/**
@@ -70,38 +67,22 @@ __attribute__((always_inline)) INLINE static void runner_iact_rt_inject(
const float r2, float *dx, const float hi, const float hj,
struct spart *restrict si, struct part *restrict pj, float a, float H) {
+ /* If the star doesn't have any neighbours, we
+ * have nothing to do here. */
+ if (si->density.wcount == 0.f) return;
+
if (si->rt_data.debug_iact_hydro_inject_prep == 0)
error(
- "Injecting energy from star that wasn't called"
- " during injection prep");
- if (pj->rt_data.debug_iact_stars_inject_prep == 0) {
+ "Injecting energy from star that wasn't called during injection prep");
- const float hig2 = hi * hi * kernel_gamma2;
- const float res = sqrtf(r2 / hig2);
- error(
- "Injecting energy into part that wasn't called"
- " during injection prep: sID %lld pID %lld r/H_s %.6f",
- si->id, pj->id, res);
- }
+ if (!si->rt_data.debug_emission_rate_set)
+ error("Injecting energy from star without setting emission rate");
si->rt_data.debug_iact_hydro_inject += 1;
si->rt_data.debug_radiation_emitted_tot += 1ULL;
pj->rt_data.debug_iact_stars_inject += 1;
pj->rt_data.debug_radiation_absorbed_tot += 1ULL;
-
- /* Attempt to catch race condition/dependency error */
- if (si->rt_data.debug_iact_hydro_inject_prep <
- si->rt_data.debug_iact_hydro_inject)
- error(
- "Star interacts with more particles during"
- " injection than during injection prep");
-
- if (pj->rt_data.debug_iact_stars_inject_prep <
- pj->rt_data.debug_iact_stars_inject)
- error(
- "Part interacts with more stars during"
- " injection than during injection prep");
}
/**
diff --git a/src/rt/debug/rt_io.h b/src/rt/debug/rt_io.h
index c9495bba1ebb5a9f7b794e91e6c37cda200d495b..8ac58f13fb7747ae9edb04dddce66f701b80800b 100644
--- a/src/rt/debug/rt_io.h
+++ b/src/rt/debug/rt_io.h
@@ -95,13 +95,8 @@ INLINE static int rt_write_particles(const struct part* parts,
"RTDebugRadAbsorbedTot", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0, parts,
rt_data.debug_radiation_absorbed_tot,
"Radiation absorbed by this part during its lifetime");
- list[7] = io_make_output_field("RTDebugStarsInjectPrepTotCounts", ULONGLONG,
- 1, UNIT_CONV_NO_UNITS, 0, parts,
- rt_data.debug_iact_stars_inject_prep_tot,
- "Total interactions with stars during "
- "injection prep during its lifetime");
- return 8;
+ return 7;
}
/**
@@ -118,8 +113,8 @@ INLINE static int rt_write_stars(const struct spart* sparts,
list[0] = io_make_output_field("RTDebugHydroIact", INT, 1, UNIT_CONV_NO_UNITS,
0, sparts, rt_data.debug_iact_hydro_inject,
- "number of interactions between this hydro "
- "particle and any star particle");
+ "number of interactions between this star "
+ "particle and any particle during injection");
list[1] = io_make_output_field(
"RTDebugEmissionRateSet", INT, 1, UNIT_CONV_NO_UNITS, 0, sparts,
rt_data.debug_emission_rate_set, "Stellar photon emission rates set?");
@@ -127,13 +122,7 @@ INLINE static int rt_write_stars(const struct spart* sparts,
"RTDebugRadEmittedTot", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0, sparts,
rt_data.debug_radiation_emitted_tot,
"Total radiation emitted during the lifetime of this star");
- list[3] = io_make_output_field("RTDebugHydroInjectPrepCountsTot", ULONGLONG,
- 1, UNIT_CONV_NO_UNITS, 0, sparts,
- rt_data.debug_iact_hydro_inject_prep_tot,
- "Total interactions with particles during "
- "injection prep during its lifetime");
-
- return 4;
+ return 3;
}
/**
@@ -153,12 +142,7 @@ INLINE static void rt_write_flavour(hid_t h_grp, hid_t h_grp_columns,
const struct rt_props* rtp) {
#if defined(HAVE_HDF5)
- if (rtp->hydro_controlled_injection) {
- io_write_attribute_s(h_grp, "RT Scheme",
- RT_IMPLEMENTATION ", hydro controlled injection");
- } else {
- io_write_attribute_s(h_grp, "RT Scheme", RT_IMPLEMENTATION);
- }
+ io_write_attribute_s(h_grp, "RT Scheme", RT_IMPLEMENTATION);
#endif
}
diff --git a/src/rt/debug/rt_properties.h b/src/rt/debug/rt_properties.h
index f26ee4bdfe96252279fefc3c770f6143bf3e549c..1c5c0b6ff1fb141f2b51a204b79b855b1d1d6c53 100644
--- a/src/rt/debug/rt_properties.h
+++ b/src/rt/debug/rt_properties.h
@@ -28,22 +28,11 @@
* @brief Properties of the debug radiative transfer model
*/
struct rt_props {
- /* Are we running with hydro or star controlled injection?
- * This is added to avoid #ifdef macros as far as possible */
- int hydro_controlled_injection;
-
- /* Do we need to run a conversion after the zeroth
- * step, but before the first step? */
- int convert_stars_after_zeroth_step;
- int convert_parts_after_zeroth_step;
-
- /* Do extended tests where we assume that all parts
- * have spart neighbours? */
- int debug_do_all_parts_have_stars_checks;
/* radiation emitted by stars this step. This is not really a property,
- * but a placeholder to sum up a global variable */
- int debug_radiation_emitted_this_step;
+ * but a placeholder to sum up a global variable. It's being reset
+ * every timestep. */
+ unsigned long long debug_radiation_emitted_this_step;
/* total radiation emitted by stars. This is not really a property,
* but a placeholder to sum up a global variable */
@@ -51,27 +40,11 @@ struct rt_props {
/* radiation absorbed by gas this step. This is not really a property,
* but a placeholder to sum up a global variable */
- int debug_radiation_absorbed_this_step;
+ unsigned long long debug_radiation_absorbed_this_step;
/* total radiation absorbed by gas. This is not really a property,
* but a placeholder to sum up a global variable */
unsigned long long debug_radiation_absorbed_tot;
-
- /* Interactions of a star with gas during injection prep this step. This is
- * not really a property, but a placeholder to sum up a global variable */
- int debug_star_injection_prep_iacts_with_parts_this_step;
-
- /* Interactions of a star with gas during injection prep. This is not
- * really a property, but a placeholder to sum up a global variable */
- unsigned long long debug_star_injection_prep_iacts_with_parts_tot;
-
- /* Interactions of a star with gas during injection prep this step. This is
- * not really a property, but a placeholder to sum up a global variable */
- int debug_part_injection_prep_iacts_with_stars_this_step;
-
- /* Interactions of a star with gas during injection prep. This is not
- * really a property, but a placeholder to sum up a global variable */
- unsigned long long debug_part_injection_prep_iacts_with_stars_tot;
};
/**
@@ -86,12 +59,6 @@ __attribute__((always_inline)) INLINE static void rt_props_print(
if (engine_rank != 0) return;
message("Radiative transfer scheme: '%s'", RT_IMPLEMENTATION);
-
- /* Print the RT properties */
- if (rtp->debug_do_all_parts_have_stars_checks)
- message("Doing extra checks assuming all parts have spart neighbours");
- else
- message("Skipping extra checks assuming all parts have spart neighbours");
}
/**
@@ -108,23 +75,11 @@ __attribute__((always_inline)) INLINE static void rt_props_init(
const struct unit_system* us, struct swift_params* params,
struct cosmology* cosmo) {
- rtp->debug_do_all_parts_have_stars_checks =
- parser_get_opt_param_int(params, "DebugRT:all_parts_have_stars", 0);
-
-#ifdef RT_HYDRO_CONTROLLED_INJECTION
- rtp->hydro_controlled_injection = 1;
-#else
- rtp->hydro_controlled_injection = 0;
-#endif
-
- /* Make sure we reset debugging counters correctly. */
- rtp->convert_parts_after_zeroth_step = 1;
- rtp->convert_stars_after_zeroth_step = 1;
-
rtp->debug_radiation_emitted_tot = 0ULL;
+ rtp->debug_radiation_emitted_this_step = 0ULL;
+
rtp->debug_radiation_absorbed_tot = 0ULL;
- rtp->debug_star_injection_prep_iacts_with_parts_tot = 0LL;
- rtp->debug_part_injection_prep_iacts_with_stars_tot = 0LL;
+ rtp->debug_radiation_absorbed_this_step = 0ULL;
/* After initialisation, print params to screen */
rt_props_print(rtp);
diff --git a/src/rt/debug/rt_stellar_emission_rate.h b/src/rt/debug/rt_stellar_emission_rate.h
deleted file mode 100644
index 348ae6435ab02a307e75a281d5d8cbd29531acff..0000000000000000000000000000000000000000
--- a/src/rt/debug/rt_stellar_emission_rate.h
+++ /dev/null
@@ -1,48 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2021 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#ifndef SWIFT_RT_STELLAR_EMISSION_RATE_DEBUG_H
-#define SWIFT_RT_STELLAR_EMISSION_RATE_DEBUG_H
-
-/**
- * @file src/rt/debug/rt_stellar_emission_rate.h
- * @brief Main header file for the debug radiative transfer scheme
- * stellar radiation emission rates related functions.
- */
-
-/**
- * @brief Main function for getting the stellar emission rate and updating it
- * in the spart.
- *
- * @param sp Star particle to work on.
- * @param age_beg Age of the stars at the beginning of the step
- * @param age_end Age of the stars at the end of the step
- * @param rt_props RT properties struct
- * @param phys_const struct holding physical constants
- * @param internal_units units struct containing internal units
- */
-
-__attribute__((always_inline)) INLINE static void rt_set_stellar_emission_rate(
- struct spart* restrict sp, double age_beg, double age_end,
- const struct rt_props* rt_props, const struct phys_const* phys_const,
- const struct unit_system* internal_units) {
-
- sp->rt_data.debug_emission_rate_set += 1;
-}
-
-#endif /* SWIFT_RT_STELLAR_EMISSION_RATE_DEBUG_H */
diff --git a/src/rt/debug/rt_struct.h b/src/rt/debug/rt_struct.h
index 03e6b9d09e5b8c623fef984441727c987cebb596..2ab50cc7513979b34a901dcf0ccfb725a295c920 100644
--- a/src/rt/debug/rt_struct.h
+++ b/src/rt/debug/rt_struct.h
@@ -32,25 +32,13 @@ struct rt_part_data {
/*! how much radiation this part received from stars during total lifetime */
unsigned long long debug_radiation_absorbed_tot;
- /*! how many interactions this part had with stars in injection prep over
- * total lifetime */
- unsigned long long debug_iact_stars_inject_prep_tot;
-
/* data to store during one time step */
- /*! how many stars this part interacted with during preparation*/
- /* Note: It's useless to write this in outputs, as it gets reset
- * at the end of every step. */
- int debug_iact_stars_inject_prep;
-
/*! how many stars this part interacted with during injection*/
/* Note: It's useless to write this in outputs, as it gets reset
* at the end of every step. */
int debug_iact_stars_inject;
- /*! called in a self/rt_injection task? */
- int debug_injection_check;
-
/*! calls from gradient interaction loop in actual function */
int debug_calls_iact_gradient_interaction;
@@ -83,10 +71,6 @@ struct rt_spart_data {
/*! how much radiation this star emitted during total lifetime */
unsigned long long debug_radiation_emitted_tot;
- /*! how many interactions this star had with parts during
- * injection prep over total lifetime */
- unsigned long long debug_iact_hydro_inject_prep_tot;
-
/* data to store during one time step */
/*! how many hydro particles this particle interacted with
@@ -99,9 +83,6 @@ struct rt_spart_data {
/*! stellar photon emisison rate computed? */
int debug_emission_rate_set;
-
- /*! called in a self/rt_injection task? */
- int debug_injection_check;
};
#endif /* SWIFT_RT_STRUCT_DEBUG_H */
diff --git a/src/rt/debug/rt_thermochemistry.h b/src/rt/debug/rt_thermochemistry.h
deleted file mode 100644
index 718926dde0bbf9fdc924bfe6ff3beae7c70bb2ff..0000000000000000000000000000000000000000
--- a/src/rt/debug/rt_thermochemistry.h
+++ /dev/null
@@ -1,49 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2020 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#ifndef SWIFT_RT_THERMOCHEMISTRY_DEBUG_H
-#define SWIFT_RT_THERMOCHEMISTRY_DEBUG_H
-
-/**
- * @file src/rt/debug/rt_thermochemistry.h
- * @brief Main header file for the debug radiative transfer scheme
- * thermochemistry related functions.
- */
-
-/**
- * @brief Main function for the thermochemistry step.
- *
- * @param p Particle to work on.
- */
-__attribute__((always_inline)) INLINE static void rt_do_thermochemistry(
- struct part *restrict p) {
-
- if (p->rt_data.debug_kicked != 1)
- error("Trying to do thermochemistry on unkicked particle %lld (count=%d)",
- p->id, p->rt_data.debug_kicked);
- if (!p->rt_data.debug_injection_done)
- error("Trying to do thermochemistry when injection step hasn't been done");
- if (!p->rt_data.debug_gradients_done)
- error("Trying to do thermochemistry when gradient step hasn't been done");
- if (!p->rt_data.debug_transport_done)
- error("Trying to do thermochemistry when transport step hasn't been done");
-
- p->rt_data.debug_thermochem_done += 1;
-}
-
-#endif /* SWIFT_RT_THERMOCHEMISTRY_DEBUG_H */
diff --git a/src/rt/none/rt.h b/src/rt/none/rt.h
index 41fafeb1f085f04576340b9bac9bfdf9fdfa47aa..6a0702813fe85ba5e0fcfc0f048f1514b456d993 100644
--- a/src/rt/none/rt.h
+++ b/src/rt/none/rt.h
@@ -80,7 +80,7 @@ __attribute__((always_inline)) INLINE static void rt_first_init_part(
/**
* @brief Initialises particle quantities that can't be set
* otherwise before the zeroth step is finished. E.g. because
- * they require the particle density to be known.
+ * they require the particle density and time step to be known.
*
* @param p particle to work on
* @param rt_props RT properties struct
@@ -120,7 +120,7 @@ __attribute__((always_inline)) INLINE static void rt_first_init_spart(
/**
* @brief Initialises particle quantities that can't be set
* otherwise before the zeroth step is finished. E.g. because
- * they require the star density to be known.
+ * they require the star density and time step to be known.
* @param sp star particle to work on
* @param time current system time
* @param star_age age of the star *at the end of the step*
@@ -306,8 +306,9 @@ __attribute__((always_inline)) INLINE static void rt_prepare_force(
* @brief Clean the allocated memory inside the RT properties struct.
*
* @param props the #rt_props.
+ * @param restart did we restart?
*/
__attribute__((always_inline)) INLINE static void rt_clean(
- struct rt_props* props) {}
+ struct rt_props* props, int restart) {}
#endif /* SWIFT_RT_NONE_H */
diff --git a/src/rt/none/rt_properties.h b/src/rt/none/rt_properties.h
index 02826ce2b992c3dfe5cae48fac4a95645e2ebc13..14b8463bb24fad7e278cbe95a9c2fb93524069c5 100644
--- a/src/rt/none/rt_properties.h
+++ b/src/rt/none/rt_properties.h
@@ -27,17 +27,7 @@
/**
* @brief Properties of the 'none' radiative transfer model
*/
-struct rt_props {
-
- /* Are we running with hydro or star controlled injection?
- * This is added to avoid #ifdef macros as far as possible */
- int hydro_controlled_injection;
-
- /* Do we need to run a conversion after the zeroth
- * step, but before the first step? */
- int convert_stars_after_zeroth_step;
- int convert_parts_after_zeroth_step;
-};
+struct rt_props {};
/**
* @brief Print the RT model.
@@ -67,15 +57,6 @@ __attribute__((always_inline)) INLINE static void rt_props_init(
const struct unit_system* us, struct swift_params* params,
struct cosmology* cosmo) {
-#ifdef RT_HYDRO_CONTROLLED_INJECTION
- rtp->hydro_controlled_injection = 1;
-#else
- rtp->hydro_controlled_injection = 0;
-#endif
-
- rtp->convert_parts_after_zeroth_step = 0;
- rtp->convert_stars_after_zeroth_step = 0;
-
/* After initialisation, print params to screen */
rt_props_print(rtp);
}
diff --git a/src/rt_active.h b/src/rt_active.h
deleted file mode 100644
index 43733ca98feda16b905216d8c96c9dbd5a2e7344..0000000000000000000000000000000000000000
--- a/src/rt_active.h
+++ /dev/null
@@ -1,142 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2021 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#ifndef SWIFT_RT_ACTIVE_H
-#define SWIFT_RT_ACTIVE_H
-
-/* Local includes. */
-#include "active.h"
-
-/**
- * @file rt_active.h
- * @brief This file contains the checks for whether radiative transfer
- * (injection) tasks should be activated for given cells, parts, or sparts.
- * The functions differ from the checks in `src/active.h` in that not only
- * time-step data is being checked, but more properties as well. So for
- * consistency, they get their own file. Finally, the functions are gathered
- * in one RT related file to concentrate all #ifdef macro shenanigans in a
- * single place as far as possible.
- */
-
-/**
- * @brief Pick whether we need to check whether an spart is active depending
- * on which version of injection we are using. This is done here to minimize
- * #ifdef macros throughout this file.
- *
- * Returns 1 if spart is active.
- *
- * @param sp star particle
- * @param e the engine
- */
-__attribute__((always_inline)) INLINE static int rt_is_spart_active_in_loop(
- struct spart *restrict sp, const struct engine *e) {
-
-#ifdef RT_HYDRO_CONTROLLED_INJECTION
- return 1; /* ignore stellar activity when gas pulls radiation from stars */
-#else
- return spart_is_active(sp, e);
-#endif
-}
-
-/**
- * @brief Pick whether we need to check whether a part is active depending
- * on which version of injection we are using. This is done here to minimize
- * #ifdef macros throughout this file.
- *
- * Returns 1 if spart is active.
- *
- * @param p the part
- * @param e the engine
- */
-__attribute__((always_inline)) INLINE static int rt_is_part_active_in_loop(
- struct part *restrict p, const struct engine *e) {
-
-#ifdef RT_HYDRO_CONTROLLED_INJECTION
- return part_is_active(p, e);
-#else
- return 1; /* ignore hydro activity when stars push radiation onto gas */
-#endif
-}
-
-/**
- * @brief Does a cell contain particles that should do RT this step?
- * This function is for a self-type interaction, where we need a cell
- * to have active hydro particles and star particles in any state.
- *
- * @param c The #cell.
- * @param e The #engine containing information about the current time.
- * @return 1 if the #cell contains at least an active particle, 0 otherwise.
- */
-__attribute__((always_inline)) INLINE static int rt_should_iact_cell(
- const struct cell *c, const struct engine *e) {
-
-#ifdef RT_HYDRO_CONTROLLED_INJECTION
- return ((cell_is_active_hydro(c, e) && (c->hydro.count > 0)) &&
- (c->stars.count > 0));
-#else
- return ((cell_is_active_stars(c, e) && (c->stars.count > 0)) &&
- (c->hydro.count > 0));
-#endif
-}
-
-/**
- * @brief Does a cell contain particles that should do RT this step?
- * This function is for a pair-type interaction, where we take stars from
- * cell ci and hydro particles from cell cj.
- *
- * @param ci First #cell.
- * @param cj Second #cell.
- * @param e The #engine containing information about the current time.
- * @return 1 if the #cell contains at least an active particle, 0 otherwise.
- */
-__attribute__((always_inline)) INLINE static int rt_should_iact_cell_pair(
- const struct cell *ci, const struct cell *cj, const struct engine *e) {
-
- if (cj == NULL) return 0;
-
-#ifdef RT_HYDRO_CONTROLLED_INJECTION
- return (cell_is_active_hydro(cj, e) && (cj->hydro.count > 0) &&
- (ci->stars.count > 0));
-#else
- return (cell_is_active_stars(ci, e) && (ci->stars.count > 0) &&
- (cj->hydro.count > 0));
-#endif
-}
-
-/**
- * @brief Do we need to unskip this cell's RT (injection) related tasks?
- * For unskipping (recursively), don't check about the star's count here:
- * Pair-type interactions don't require stars in every cell. This way, we
- * can include the check for star count > 0 there on a pair by pair basis.
- *
- * @param c The #cell.
- * @param e The #engine containing information about the current time.
- * @return 1 if the #cell needs to activate tasks
- */
-__attribute__((always_inline)) INLINE static int rt_should_do_unskip_cell(
- const struct cell *c, const struct engine *e) {
- /* whether it's hydro controlled or not, we need to check for hydro
- * activity at the top level. We also need to check for star activity
- * so we can activate the rt_in implicit tasks to catch dependencies
- * before the injection and not be overwritten by work in star density
- * ghosts. */
- return ((cell_is_active_hydro(c, e) && (c->hydro.count > 0)) ||
- cell_is_active_stars(c, e));
-}
-
-#endif /* SWIFT_RT_ACTIVE_H */
diff --git a/src/runner.h b/src/runner.h
index 6f9bbe582438704153b770dad6628a8d3283808f..af450ba402727d3fe2da3c713d481853e2121d79 100644
--- a/src/runner.h
+++ b/src/runner.h
@@ -48,7 +48,6 @@ struct task;
#define TASK_LOOP_STARS_PREP2 10
#define TASK_LOOP_RT_GRADIENT 11
#define TASK_LOOP_RT_TRANSPORT 12
-#define TASK_LOOP_RT_INJECT 13
/**
* @brief A struct representing a runner's thread and its data.
diff --git a/src/runner_doiact_functions_rt.h b/src/runner_doiact_functions_rt.h
deleted file mode 100644
index 1e73acf20b5ea4511646f1c6979f7d9037e57914..0000000000000000000000000000000000000000
--- a/src/runner_doiact_functions_rt.h
+++ /dev/null
@@ -1,809 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
- * 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- * 2020 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-
-/* Before including this file, define FUNCTION, which is the
- name of the interaction function. This creates the interaction functions
- runner_dopair_FUNCTION, runner_dopair_FUNCTION_naive, runner_doself_FUNCTION,
- and runner_dosub_FUNCTION calling the pairwise interaction function
- runner_iact_FUNCTION. */
-
-#include "rt.h"
-#include "rt_active.h"
-#include "runner_doiact_rt.h"
-
-/**
- * @brief Function for self-type interaction between stars and hydro particles
- *
- * @param r runner task
- * @param c cell
- * @param timer 1 if the time is to be recorded.
- */
-void DOSELF1_RT(struct runner *r, struct cell *c, int timer) {
-
- TIMER_TIC;
-
- const struct engine *e = r->e;
-
- /* Cosmological terms */
- const struct cosmology *cosmo = e->cosmology;
- const float a = cosmo->a;
- const float H = cosmo->H;
-
- /* Anything to do here? */
- if (c->hydro.count == 0 || c->stars.count == 0) return;
-
-#ifdef SWIFT_DEBUG_CHECKS
- /* Drift the cell to the current timestep if needed. */
- if (!cell_are_spart_drifted(c, r->e))
- error("Interacting undrifted cell (spart): %lld", c->cellID);
- if (!cell_are_part_drifted(c, r->e))
- error("Interacting undrifted cell (part): %lld", c->cellID);
-#endif
-
- struct spart *restrict sparts = c->stars.parts;
- struct part *restrict parts = c->hydro.parts;
-
- const int scount = c->stars.count;
- const int count = c->hydro.count;
-
- /* Loop over the sparts in cell */
- for (int sid = 0; sid < scount; sid++) {
-
- /* Get a hold of the ith spart in c. */
- struct spart *restrict si = &sparts[sid];
-
- /* Skip inhibited particles. */
- if (spart_is_inhibited(si, e)) continue;
-
- /* Skip inactive particles */
- if (!rt_is_spart_active_in_loop(si, e)) continue;
-
- const float hi = si->h;
- const float hig2 = hi * hi * kernel_gamma2;
- const float six[3] = {(float)(si->x[0] - c->loc[0]),
- (float)(si->x[1] - c->loc[1]),
- (float)(si->x[2] - c->loc[2])};
-
- /* Loop over the parts in cell */
- for (int pid = 0; pid < count; pid++) {
-
- /* Get a pointer to the jth particle. */
- struct part *restrict pj = &parts[pid];
- const float hj = pj->h;
-
- /* Skip inhibited particles. */
- if (part_is_inhibited(pj, e)) continue;
-
- /* Skip inactive particles. */
- if (!rt_is_part_active_in_loop(pj, e)) continue;
-
- /* Compute the pairwise distance. */
- const float pjx[3] = {(float)(pj->x[0] - c->loc[0]),
- (float)(pj->x[1] - c->loc[1]),
- (float)(pj->x[2] - c->loc[2])};
- float dx[3] = {six[0] - pjx[0], six[1] - pjx[1], six[2] - pjx[2]};
- const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
-
-#ifdef SWIFT_DEBUG_CHECKS
- /* Check that particles have been drifted to the current time */
- if (pj->ti_drift != e->ti_current)
- error("Particle pj not drifted to current time");
-
- if (r2 < hig2 && e->rt_props->hydro_controlled_injection) {
-
- const integertime_t ti_current = e->ti_current;
- const integertime_t pti_end =
- get_integer_time_end(ti_current, pj->time_bin);
- const integertime_t sti_end =
- get_integer_time_end(ti_current, si->time_bin);
-
- if (sti_end < ti_current)
- error(
- "s-particle in an impossible time-zone! sp->ti_end=%lld "
- "e->ti_current=%lld",
- sti_end, ti_current);
- if (pti_end < ti_current)
- error(
- "particle in an impossible time-zone! p->ti_end=%lld "
- "e->ti_current=%lld",
- pti_end, ti_current);
- if (pti_end > sti_end)
- message(
- "WARNING: Got star that whose time step ends before the "
- "interacting "
- "hydro particle's time step ends. This needs to be dealt with.");
- }
-#endif
-
- if (r2 < hig2) {
- IACT_RT(r2, dx, hi, hj, si, pj, a, H);
- }
- }
- }
-
- if (timer) TIMER_TOC(TIMER_DOSELF_RT);
-}
-
-/**
- * @brief Function for non-symmetric pair-type interaction between stars
- * and hydro particles. Will interact star particles of cell i
- * with hydro particles of cell j.
- *
- *
- * @param r runner task
- * @param ci the first cell, where we take star particles from
- * @param cj the second cell, where we take hydro particles from
- */
-void DOPAIR1_NONSYM_RT_NAIVE(struct runner *r, struct cell *ci,
- struct cell *cj) {
-
- TIMER_TIC;
-
- const struct engine *e = r->e;
-
- /* Cosmological terms */
- const struct cosmology *cosmo = e->cosmology;
- const float a = cosmo->a;
- const float H = cosmo->H;
-
- const int scount_i = ci->stars.count;
- const int count_j = cj->hydro.count;
- struct spart *restrict sparts_i = ci->stars.parts;
- struct part *restrict parts_j = cj->hydro.parts;
-
- /* Get the relative distance between the pairs, wrapping. */
- double shift[3] = {0.0, 0.0, 0.0};
- for (int k = 0; k < 3; k++) {
- if (cj->loc[k] - ci->loc[k] < -e->s->dim[k] / 2)
- shift[k] = e->s->dim[k];
- else if (cj->loc[k] - ci->loc[k] > e->s->dim[k] / 2)
- shift[k] = -e->s->dim[k];
- }
-
- /* Loop over the sparts in ci. */
- for (int sid = 0; sid < scount_i; sid++) {
-
- /* Get a hold of the ith spart in ci. */
- struct spart *restrict si = &sparts_i[sid];
-
- /* Skip inhibited particles. */
- if (spart_is_inhibited(si, e)) continue;
-
- /* Skip inactive particles */
- if (!rt_is_spart_active_in_loop(si, e)) continue;
-
- const float hi = si->h;
- const float hig2 = hi * hi * kernel_gamma2;
- const float six[3] = {(float)(si->x[0] - (cj->loc[0] + shift[0])),
- (float)(si->x[1] - (cj->loc[1] + shift[1])),
- (float)(si->x[2] - (cj->loc[2] + shift[2]))};
-
- /* Loop over the parts in cj. */
- for (int pjd = 0; pjd < count_j; pjd++) {
-
- /* Get a pointer to the jth particle. */
- struct part *restrict pj = &parts_j[pjd];
- const float hj = pj->h;
-
- /* Skip inhibited particles. */
- if (part_is_inhibited(pj, e)) continue;
-
- /* Skip inactive particles. */
- if (!rt_is_part_active_in_loop(pj, e)) continue;
-
- /* Compute the pairwise distance. */
- const float pjx[3] = {(float)(pj->x[0] - cj->loc[0]),
- (float)(pj->x[1] - cj->loc[1]),
- (float)(pj->x[2] - cj->loc[2])};
- float dx[3] = {six[0] - pjx[0], six[1] - pjx[1], six[2] - pjx[2]};
- const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
-
-#ifdef RT_DEBUG
- if (r2 < hig2 && e->rt_props->hydro_controlled_injection) {
-
- const integertime_t ti_current = e->ti_current;
- const integertime_t pti_end =
- get_integer_time_end(ti_current, pj->time_bin);
- const integertime_t sti_end =
- get_integer_time_end(ti_current, si->time_bin);
-
- if (sti_end < ti_current)
- error(
- "s-particle in an impossible time-zone! sp->ti_end=%lld "
- "e->ti_current=%lld",
- sti_end, ti_current);
- if (pti_end < ti_current)
- error(
- "particle in an impossible time-zone! p->ti_end=%lld "
- "e->ti_current=%lld",
- pti_end, ti_current);
- if (pti_end > sti_end)
- message(
- "WARNING: Got star that whose time step ends before the "
- "interacting "
- "hydro particle's time step ends. This needs to be dealt with.");
- }
-#endif
-
- if (r2 < hig2) {
- IACT_RT(r2, dx, hi, hj, si, pj, a, H);
- }
-
- } /* loop over the parts in cj. */
- } /* loop over the parts in ci. */
-}
-
-/**
- * @brief Function for pair-type interaction between stars
- * and hydro particles. Will interact hydro particles of cell i
- * with star particles of cell j.
- *
- * @param r #runner task
- * @param ci the first #cell
- * @param cj the second #cell
- * @param timer 1 if the time is to be recorded.
- */
-void DO_SYM_PAIR1_RT(struct runner *r, struct cell *ci, struct cell *cj,
- const int sid, const double *shift) {
-
- TIMER_TIC;
-
- const struct engine *e = r->e;
-
- /* Cosmological terms */
- const struct cosmology *cosmo = e->cosmology;
- const float a = cosmo->a;
- const float H = cosmo->H;
-
- /* Get the cutoff shift. */
- double rshift = 0.0;
- for (int k = 0; k < 3; k++) rshift += shift[k] * runner_shift[sid][k];
-
- const int do_ci_stars =
- (ci->nodeID == e->nodeID) && rt_should_iact_cell_pair(ci, cj, e);
- const int do_cj_stars =
- (cj->nodeID == e->nodeID) && rt_should_iact_cell_pair(cj, ci, e);
-
- if (do_ci_stars) {
-
- /* Pick-out the sorted lists. */
- const struct sort_entry *restrict sort_j = cell_get_hydro_sorts(cj, sid);
- const struct sort_entry *restrict sort_i = cell_get_stars_sorts(ci, sid);
-
- /* Get some other useful values. */
- const double hi_max = ci->stars.h_max * kernel_gamma - rshift;
- const int count_i = ci->stars.count;
- const int count_j = cj->hydro.count;
- struct spart *restrict sparts_i = ci->stars.parts;
- struct part *restrict parts_j = cj->hydro.parts;
- const double dj_min = sort_j[0].d;
- const float dx_max = (ci->stars.dx_max_sort + cj->hydro.dx_max_sort);
- const float hydro_dx_max_rshift = cj->hydro.dx_max_sort - rshift;
-
- /* TODO: maybe change the order of the loops for better performance? */
- /* Loop over the sparts in ci. */
- for (int pid = count_i - 1;
- pid >= 0 && sort_i[pid].d + hi_max + dx_max > dj_min; pid--) {
-
- /* Get a hold of the ith spart in ci. */
- struct spart *restrict spi = &sparts_i[sort_i[pid].i];
- const float hi = spi->h;
-
- /* Skip inhibited particles */
- if (spart_is_inhibited(spi, e)) continue;
-
- /* Skip inactive particles */
- if (!rt_is_spart_active_in_loop(spi, e)) continue;
-
- /* Compute distance from the other cell. */
- const double px[3] = {spi->x[0], spi->x[1], spi->x[2]};
- float dist = px[0] * runner_shift[sid][0] + px[1] * runner_shift[sid][1] +
- px[2] * runner_shift[sid][2];
-
- /* Is there anything we need to interact with ? */
- const double di = dist + hi * kernel_gamma + hydro_dx_max_rshift;
- if (di < dj_min) continue;
-
- /* Get some additional information about pi */
- const float hig2 = hi * hi * kernel_gamma2;
- const float pix = spi->x[0] - (cj->loc[0] + shift[0]);
- const float piy = spi->x[1] - (cj->loc[1] + shift[1]);
- const float piz = spi->x[2] - (cj->loc[2] + shift[2]);
-
- /* Loop over the parts in cj. */
- for (int pjd = 0; pjd < count_j && sort_j[pjd].d < di; pjd++) {
-
- /* Recover pj */
- struct part *pj = &parts_j[sort_j[pjd].i];
-
- /* Skip inhibited particles. */
- if (part_is_inhibited(pj, e)) continue;
-
- /* Skip inactive particles. */
- if (!rt_is_part_active_in_loop(pj, e)) continue;
-
- const float hj = pj->h;
- const float pjx = pj->x[0] - cj->loc[0];
- const float pjy = pj->x[1] - cj->loc[1];
- const float pjz = pj->x[2] - cj->loc[2];
-
- /* Compute the pairwise distance. */
- float dx[3] = {pix - pjx, piy - pjy, piz - pjz};
- const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
-
-#ifdef SWIFT_DEBUG_CHECKS
- /* Check that particles have been drifted to the current time */
- if (spi->ti_drift != e->ti_current)
- error("Particle spi not drifted to current time");
- if (pj->ti_drift != e->ti_current)
- error("Particle pj not drifted to current time");
-
- if (r2 < hig2 && e->rt_props->hydro_controlled_injection) {
-
- const integertime_t ti_current = e->ti_current;
- const integertime_t pti_end =
- get_integer_time_end(ti_current, pj->time_bin);
- const integertime_t sti_end =
- get_integer_time_end(ti_current, spi->time_bin);
-
- if (sti_end < ti_current)
- error(
- "s-particle in an impossible time-zone! sp->ti_end=%lld "
- "e->ti_current=%lld",
- sti_end, ti_current);
- if (pti_end < ti_current)
- error(
- "particle in an impossible time-zone! p->ti_end=%lld "
- "e->ti_current=%lld",
- pti_end, ti_current);
- if (pti_end > sti_end)
- message(
- "WARNING: Got star that whose time step ends before the "
- "interacting "
- "hydro particle's time step ends. This needs to be dealt "
- "with.");
- }
-
-#endif
-
- /* Hit or miss? */
- if (r2 < hig2) {
- IACT_RT(r2, dx, hi, hj, spi, pj, a, H);
- }
- } /* loop over the parts in cj. */
- } /* loop over the parts in ci. */
- } /* do_ci_stars */
-
- if (do_cj_stars) {
- /* Pick-out the sorted lists. */
- const struct sort_entry *restrict sort_i = cell_get_hydro_sorts(ci, sid);
- const struct sort_entry *restrict sort_j = cell_get_stars_sorts(cj, sid);
-
- /* Get some other useful values. */
- const double hj_max = cj->stars.h_max * kernel_gamma;
- const int count_i = ci->hydro.count;
- const int count_j = cj->stars.count;
- struct part *restrict parts_i = ci->hydro.parts;
- struct spart *restrict sparts_j = cj->stars.parts;
- const double di_max = sort_i[count_i - 1].d - rshift;
- const float dx_max = (ci->hydro.dx_max_sort + cj->stars.dx_max_sort);
- const float hydro_dx_max_rshift = ci->hydro.dx_max_sort - rshift;
-
- /* TODO: maybe change the order of the loops for better performance? */
- /* Loop over the sparts in cj. */
- for (int pjd = 0; pjd < count_j && sort_j[pjd].d - hj_max - dx_max < di_max;
- pjd++) {
-
- /* Get a hold of the jth spart in cj. */
- struct spart *spj = &sparts_j[sort_j[pjd].i];
- const float hj = spj->h;
-
- /* Skip inhibited particles */
- if (spart_is_inhibited(spj, e)) continue;
-
- /* Skip inactive particles */
- if (!rt_is_spart_active_in_loop(spj, e)) continue;
-
- /* Compute distance from the other cell. */
- const double px[3] = {spj->x[0], spj->x[1], spj->x[2]};
- float dist = px[0] * runner_shift[sid][0] + px[1] * runner_shift[sid][1] +
- px[2] * runner_shift[sid][2];
-
- /* Is there anything we need to interact with ? */
- const double dj = dist - hj * kernel_gamma - hydro_dx_max_rshift;
- if (dj - rshift > di_max) continue;
-
- /* Get some additional information about pj */
- const float hjg2 = hj * hj * kernel_gamma2;
- const float pjx = spj->x[0] - cj->loc[0];
- const float pjy = spj->x[1] - cj->loc[1];
- const float pjz = spj->x[2] - cj->loc[2];
-
- /* Loop over the parts in ci. */
- for (int pid = count_i - 1; pid >= 0 && sort_i[pid].d > dj; pid--) {
-
- /* Recover pi */
- struct part *pi = &parts_i[sort_i[pid].i];
-
- /* Skip inhibited particles. */
- if (part_is_inhibited(pi, e)) continue;
-
- /* Skip inactive particles. */
- if (!rt_is_part_active_in_loop(pi, e)) continue;
-
- const float hi = pi->h;
- const float pix = pi->x[0] - (cj->loc[0] + shift[0]);
- const float piy = pi->x[1] - (cj->loc[1] + shift[1]);
- const float piz = pi->x[2] - (cj->loc[2] + shift[2]);
-
- /* Compute the pairwise distance. */
- float dx[3] = {pjx - pix, pjy - piy, pjz - piz};
- const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
-
-#ifdef SWIFT_DEBUG_CHECKS
- /* Check that particles have been drifted to the current time */
- if (pi->ti_drift != e->ti_current)
- error("Particle pi not drifted to current time");
- if (spj->ti_drift != e->ti_current)
- error("Particle spj not drifted to current time");
-
- if (r2 < hjg2 && e->rt_props->hydro_controlled_injection) {
-
- const integertime_t ti_current = e->ti_current;
- const integertime_t pti_end =
- get_integer_time_end(ti_current, pi->time_bin);
- const integertime_t sti_end =
- get_integer_time_end(ti_current, spj->time_bin);
-
- if (sti_end < ti_current)
- error(
- "s-particle in an impossible time-zone! sp->ti_end=%lld "
- "e->ti_current=%lld",
- sti_end, ti_current);
- if (pti_end < ti_current)
- error(
- "particle in an impossible time-zone! p->ti_end=%lld "
- "e->ti_current=%lld",
- pti_end, ti_current);
- if (pti_end > sti_end)
- message(
- "WARNING: Got star that whose time step ends before the "
- "interacting "
- "hydro particle's time step ends. This needs to be dealt "
- "with.");
- }
-#endif
-
- /* Hit or miss? */
- if (r2 < hjg2) {
- IACT_RT(r2, dx, hj, hi, spj, pi, a, H);
- }
- } /* loop over the parts in ci. */
- } /* loop over the parts in cj. */
- } /* Cell cj is active */
-
- TIMER_TOC(TIMER_DOPAIR_RT);
-}
-
-/**
- * @brief Function for pair-type interaction between stars
- * and hydro particles. Will interact hydro particles of cell i
- * with star particles of cell j.
- *
- * @param r #runner task
- * @param ci the first #cell
- * @param cj the second #cell
- * @param timer 1 if the time is to be recorded.
- */
-void DOPAIR1_RT_NAIVE(struct runner *r, struct cell *ci, struct cell *cj,
- int timer) {
-
- TIMER_TIC;
- const struct engine *restrict e = r->e;
-
- const int do_stars_in_ci =
- (cj->nodeID == r->e->nodeID) && rt_should_iact_cell_pair(ci, cj, e);
- if (do_stars_in_ci) DOPAIR1_NONSYM_RT_NAIVE(r, ci, cj);
-
- const int do_stars_in_cj =
- (ci->nodeID == r->e->nodeID) && rt_should_iact_cell_pair(cj, ci, e);
- if (do_stars_in_cj) DOPAIR1_NONSYM_RT_NAIVE(r, cj, ci);
-
- if (timer) TIMER_TOC(TIMER_DOPAIR_RT);
-}
-
-/**
- * @brief Determine which version of DOSELF1_RT needs to be called
- *
- * @param r #runner
- * @param c #cell c
- * @param timer 1 if the time is to be recorded.
- */
-void DOSELF1_BRANCH_RT(struct runner *r, struct cell *c, int timer) {
-
-#ifdef RT_DEBUG
- /* Before an early exit, loop over all parts and sparts in this cell
- * and mark that we checked these particles */
-
- const struct engine *e = r->e;
-
- if (c->hydro.count > 0) {
- struct part *restrict parts = c->hydro.parts;
-
- /* Loop over the parts in cell */
- for (int pid = 0; pid < c->hydro.count; pid++) {
-
- /* Get a pointer to the jth particle. */
- struct part *restrict pj = &parts[pid];
-
- /* Skip inhibited particles. */
- if (part_is_inhibited(pj, e)) continue;
-
- /* Skip inactive particles. */
- if (!rt_is_part_active_in_loop(pj, e)) continue;
-
- rt_debugging_check_injection_part(pj, e->rt_props);
- }
- }
-
- if (c->stars.count > 0) {
- struct spart *restrict sparts = c->stars.parts;
-
- /* Loop over the parts in cell */
- for (int sid = 0; sid < c->stars.count; sid++) {
-
- /* Get a pointer to the ith spart. */
- struct spart *restrict si = &sparts[sid];
-
- /* Skip inhibited particles. */
- if (spart_is_inhibited(si, e)) continue;
-
- /* Skip inactive particles */
- if (!rt_is_spart_active_in_loop(si, e)) continue;
-
- rt_debugging_check_injection_spart(si, e->rt_props);
- }
- }
-#endif
-
- /* early exit? */
- if (!rt_should_iact_cell(c, r->e)) return;
- DOSELF1_RT(r, c, timer);
-}
-
-/**
- * @brief Determine which version of DOPAIR1_RT needs to be called
- *
- * @param r #runner
- * @param ci The first #cell
- * @param cj The second #cell
- * @param timer 1 if the time is to be recorded.
- */
-void DOPAIR1_BRANCH_RT(struct runner *r, struct cell *ci, struct cell *cj,
- int timer) {
-
- const struct engine *restrict e = r->e;
-
- /* Get the sort ID. */
- double shift[3] = {0.0, 0.0, 0.0};
- const int sid = space_getsid(e->s, &ci, &cj, shift);
-
- const int do_stars_ci =
- (cj->nodeID == r->e->nodeID) && rt_should_iact_cell_pair(ci, cj, e);
- const int do_stars_cj =
- (ci->nodeID == r->e->nodeID) && rt_should_iact_cell_pair(cj, ci, e);
-
- /* Anything to do here? */
- if (!do_stars_ci && !do_stars_cj) return;
-
-#ifdef SWIFT_DEBUG_CHECKS
- if (do_stars_ci) {
-
- /* Check that cells are drifted. */
- if (!cell_are_spart_drifted(ci, e))
- error("Interacting undrifted stars in cells i=%12lld, j=%12lld",
- ci->cellID, cj->cellID);
- if (!cell_are_part_drifted(cj, e))
- error("Interacting undrifted hydro in cells i=%12lld, j=%12lld",
- ci->cellID, cj->cellID);
-
- /* Have the cells been sorted? */
- if (!(ci->stars.sorted & (1 << sid)) ||
- ci->stars.dx_max_sort_old > space_maxreldx * ci->dmin)
- error("Interacting unsorted cells: %lld %d", ci->cellID, sid);
- if (!(cj->hydro.sorted & (1 << sid)) ||
- cj->hydro.dx_max_sort_old > space_maxreldx * cj->dmin)
- error("Interacting unsorted cells: %lld %lld", ci->cellID, cj->cellID);
- }
-
- if (do_stars_cj) {
-
- /* Check that cells are drifted. */
- if (!cell_are_spart_drifted(cj, e)) {
- error("Interacting undrifted stars in cells j=%12lld, i=%12lld",
- cj->cellID, ci->cellID);
- }
- if (!cell_are_part_drifted(ci, e))
- error("Interacting undrifted hydro in cells j=%12lld, i=%12lld",
- cj->cellID, ci->cellID);
-
- /* Have the cells been sorted? */
- if (!(ci->hydro.sorted & (1 << sid)) ||
- ci->hydro.dx_max_sort_old > space_maxreldx * ci->dmin)
- error("Interacting unsorted cells: %lld %d", cj->cellID, sid);
-
- if ((!(cj->stars.sorted & (1 << sid)) ||
- cj->stars.dx_max_sort_old > space_maxreldx * cj->dmin))
- error("Interacting unsorted cells: %lld %lld", ci->cellID, cj->cellID);
- }
-#endif /* SWIFT_DEBUG_CHECKS */
-
- if (do_stars_ci || do_stars_cj) {
-#ifdef SWIFT_USE_NAIVE_INTERACTIONS_RT
- DOPAIR1_RT_NAIVE(r, ci, cj, 1);
-#else
- DO_SYM_PAIR1_RT(r, ci, cj, sid, shift);
-#endif
- }
-}
-
-/**
- * @brief Compute grouped sub-cell interactions for self tasks
- *
- * @param r The #runner.
- * @param ci The first #cell.
- * @param gettimer Do we have a timer ?
- */
-void DOSUB_SELF1_RT(struct runner *r, struct cell *c, int timer) {
-
- TIMER_TIC;
-
-#ifdef SWIFT_DEBUG_CHECKS
- if (c->nodeID != engine_rank)
- error("This function should not be called on foreign cells");
-#endif
-
- /* Should we even bother? */
- if (!rt_should_iact_cell(c, r->e)) {
-
-#ifdef RT_DEBUG
- /* Before an early exit, loop over all parts and sparts in this cell
- * and mark that we checked these particles */
-
- const struct engine *e = r->e;
-
- if (c->hydro.count > 0) {
- struct part *restrict parts = c->hydro.parts;
- const int count = c->hydro.count;
-
- /* Loop over the parts in cell */
- for (int pid = 0; pid < count; pid++) {
-
- /* Get a pointer to the jth particle. */
- struct part *restrict pj = &parts[pid];
-
- /* Skip inhibited particles. */
- if (part_is_inhibited(pj, e)) continue;
-
- /* Skip inactive particles. */
- if (!rt_is_part_active_in_loop(pj, e)) continue;
-
- rt_debugging_check_injection_part(pj, e->rt_props);
- }
- }
-
- if (c->stars.count > 0) {
- struct spart *restrict sparts = c->stars.parts;
- const int scount = c->stars.count;
-
- /* Loop over the parts in cell */
- for (int sid = 0; sid < scount; sid++) {
-
- /* Get a pointer to the ith spart. */
- struct spart *restrict si = &sparts[sid];
-
- /* Skip inhibited particles. */
- if (spart_is_inhibited(si, e)) continue;
-
- /* Skip inactive particles */
- if (!rt_is_spart_active_in_loop(si, e)) continue;
-
- rt_debugging_check_injection_spart(si, e->rt_props);
- }
- }
-#endif
-
- /* exit early if there is nothing to do */
- return;
- }
-
- /* Recurse? */
- if (cell_can_recurse_in_self_stars_task(c)) {
-
- /* Loop over all progeny. */
- for (int k = 0; k < 8; k++)
- if (c->progeny[k] != NULL) {
- DOSUB_SELF1_RT(r, c->progeny[k], 0);
- for (int j = k + 1; j < 8; j++)
- if (c->progeny[j] != NULL)
- DOSUB_PAIR1_RT(r, c->progeny[k], c->progeny[j], 0);
- }
- }
-
- /* Otherwise, compute self-interaction. */
- else {
- DOSELF1_BRANCH_RT(r, c, 0);
- }
-
- if (timer) TIMER_TOC(TIMER_DOSUB_SELF_RT);
-}
-
-/**
- * @brief Compute grouped sub-cell interactions for pair tasks
- *
- * @param r The #runner.
- * @param ci The first #cell.
- * @param cj The second #cell.
- * @param gettimer Do we have a timer ?
- */
-void DOSUB_PAIR1_RT(struct runner *r, struct cell *ci, struct cell *cj,
- int timer) {
-
- TIMER_TIC;
-
- struct space *s = r->e->s;
- const struct engine *e = r->e;
-
- /* Should we even bother? */
- const int should_do_ci = rt_should_iact_cell_pair(ci, cj, e);
- const int should_do_cj = rt_should_iact_cell_pair(cj, ci, e);
- if (!should_do_ci && !should_do_cj) return;
-
- /* Get the type of pair and flip ci/cj if needed. */
- double shift[3];
- const int sid = space_getsid(s, &ci, &cj, shift);
-
- /* Recurse? */
- if (cell_can_recurse_in_pair_stars_task(ci, cj) &&
- cell_can_recurse_in_pair_stars_task(cj, ci)) {
- struct cell_split_pair *csp = &cell_split_pairs[sid];
- for (int k = 0; k < csp->count; k++) {
- const int pid = csp->pairs[k].pid;
- const int pjd = csp->pairs[k].pjd;
- if (ci->progeny[pid] != NULL && cj->progeny[pjd] != NULL)
- DOSUB_PAIR1_RT(r, ci->progeny[pid], cj->progeny[pjd], 0);
- }
- }
-
- /* Otherwise, compute the pair directly. */
- else {
-
- /* do full checks again, space_getsid() might swap ci/cj pointers */
- const int do_ci_stars =
- (cj->nodeID == e->nodeID) && rt_should_iact_cell_pair(ci, cj, e);
- const int do_cj_stars =
- (ci->nodeID == e->nodeID) && rt_should_iact_cell_pair(cj, ci, e);
-
- if (do_ci_stars || do_cj_stars) DOPAIR1_BRANCH_RT(r, ci, cj, 0);
- }
-
- if (timer) TIMER_TOC(TIMER_DOSUB_PAIR_RT);
-}
diff --git a/src/runner_doiact_functions_stars.h b/src/runner_doiact_functions_stars.h
index 88ac45b4465021587158ca71e1b52774d04d0a2a..effa7751291c01eac2db6d1345a2651201451000 100644
--- a/src/runner_doiact_functions_stars.h
+++ b/src/runner_doiact_functions_stars.h
@@ -26,6 +26,17 @@
#include "runner_doiact_stars.h"
+#ifdef RT_NONE
+#define WITH_RT 0
+#else
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) || \
+ (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
+#define WITH_RT (e->policy & engine_policy_rt)
+#else
+#define WITH_RT 0
+#endif
+#endif
+
/**
* @brief Calculate the number density of #part around the #spart
*
@@ -60,13 +71,8 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
#if (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
struct xpart *restrict xparts = c->hydro.xparts;
#endif
-#if (!(defined RT_NONE) && (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK))
- /* Don't exit early if star isn't active for feedback
- * when we're running with RT */
- const int with_rt = (e->policy & engine_policy_rt);
-#else
- const int with_rt = 0;
-#endif
+
+ const int with_rt = WITH_RT;
/* Loop over the sparts in ci. */
for (int sid = 0; sid < scount; sid++) {
@@ -131,10 +137,16 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, si, pj, xpj, cosmo,
e->hydro_properties,
e->feedback_props, ti_current);
+#endif
}
if (r2 < hig2 && with_rt) {
+ /* If we're running RT, we don't care whether star is active for
+ * feedback, just that the star is active. */
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_rt_injection_prep(r2, dx, hi, hj, si, pj, cosmo,
e->rt_props);
+#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
+ runner_iact_rt_inject(r2, dx, hi, hj, si, pj, a, H);
#endif
}
} /* loop over the parts in ci. */
@@ -180,13 +192,7 @@ void DO_NONSYM_PAIR1_STARS_NAIVE(struct runner *r, struct cell *restrict ci,
#if (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
struct xpart *restrict xparts_j = cj->hydro.xparts;
#endif
-#if (!(defined RT_NONE) && (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK))
- /* Don't exit early if star isn't active for feedback
- * when we're running with RT */
- const int with_rt = (e->policy & engine_policy_rt);
-#else
- const int with_rt = 0;
-#endif
+ const int with_rt = WITH_RT;
/* Get the relative distance between the pairs, wrapping. */
double shift[3] = {0.0, 0.0, 0.0};
@@ -261,10 +267,16 @@ void DO_NONSYM_PAIR1_STARS_NAIVE(struct runner *r, struct cell *restrict ci,
runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, si, pj, xpj, cosmo,
e->hydro_properties,
e->feedback_props, ti_current);
+#endif
}
if (r2 < hig2 && with_rt) {
+ /* If we're running RT, we don't care whether star is active for
+ * feedback, just that the star is active. */
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_rt_injection_prep(r2, dx, hi, hj, si, pj, cosmo,
e->rt_props);
+#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
+ runner_iact_rt_inject(r2, dx, hi, hj, si, pj, a, H);
#endif
}
} /* loop over the parts in cj. */
@@ -310,14 +322,7 @@ void DO_SYM_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
const int do_cj_stars = (ci->nodeID == e->nodeID) && (cj->stars.count != 0) &&
(ci->hydro.count != 0) && cell_is_active_stars(cj, e);
#endif
-
-#if (!(defined RT_NONE) && (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK))
- /* Don't exit early if star isn't active for feedback
- * when we're running with RT */
- const int with_rt = (e->policy & engine_policy_rt);
-#else
- const int with_rt = 0;
-#endif
+ const int with_rt = WITH_RT;
if (do_ci_stars) {
@@ -457,10 +462,16 @@ void DO_SYM_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, xpj, cosmo,
e->hydro_properties,
e->feedback_props, ti_current);
+#endif
}
if (r2 < hig2 && with_rt) {
+ /* If we're running RT, we don't care whether star is active for
+ * feedback, just that the star is active. */
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_rt_injection_prep(r2, dx, hi, hj, spi, pj, cosmo,
e->rt_props);
+#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
+ runner_iact_rt_inject(r2, dx, hi, hj, spi, pj, a, H);
#endif
}
} /* loop over the parts in cj. */
@@ -605,10 +616,16 @@ void DO_SYM_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
runner_iact_nonsym_feedback_apply(r2, dx, hj, hi, spj, pi, xpi, cosmo,
e->hydro_properties,
e->feedback_props, ti_current);
+#endif
}
if (r2 < hjg2 && with_rt) {
+ /* If we're running RT, we don't care whether star is active for
+ * feedback, just that the star is active. */
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_rt_injection_prep(r2, dx, hj, hi, spj, pi, cosmo,
e->rt_props);
+#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
+ runner_iact_rt_inject(r2, dx, hj, hi, spj, pi, a, H);
#endif
}
} /* loop over the parts in ci. */
@@ -729,6 +746,8 @@ void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, NULL,
cosmo, e->feedback_props,
e->ti_current);
+ runner_iact_nonsym_rt_injection_prep(r2, dx, hi, hj, spi, pj, cosmo,
+ e->rt_props);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
error("No subset feedback iact functions do (or should) exist!");
/* runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, xpj,
@@ -790,6 +809,8 @@ void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, NULL,
cosmo, e->feedback_props,
e->ti_current);
+ runner_iact_nonsym_rt_injection_prep(r2, dx, hi, hj, spi, pj, cosmo,
+ e->rt_props);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
error("No subset feedback iact functions do (or should) exist!");
/* runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, xpj,
@@ -886,6 +907,8 @@ void DOPAIR1_SUBSET_STARS_NAIVE(struct runner *r, struct cell *restrict ci,
runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, NULL,
cosmo, e->feedback_props,
e->ti_current);
+ runner_iact_nonsym_rt_injection_prep(r2, dx, hi, hj, spi, pj, cosmo,
+ e->rt_props);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
error("No subset feedback iact functions do (or should) exist! .");
/* runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, xpj,
@@ -971,6 +994,8 @@ void DOSELF1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
runner_iact_nonsym_feedback_density(r2, dx, hi, pj->h, spi, pj, NULL,
cosmo, e->feedback_props,
e->ti_current);
+ runner_iact_nonsym_rt_injection_prep(r2, dx, hi, pj->h, spi, pj, cosmo,
+ e->rt_props);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
error("No subset feedback iact functions do (or should) exist!");
/* runner_iact_nonsym_feedback_apply(r2, dx, hi, pj->h, spi, pj, xpj, */
@@ -1427,3 +1452,5 @@ void DOSUB_SELF1_STARS(struct runner *r, struct cell *ci, int gettimer) {
TIMER_TOC(TIMER_DOSUB_SELF_STARS);
}
+
+#undef WITH_RT
diff --git a/src/runner_doiact_rt.c b/src/runner_doiact_rt.c
deleted file mode 100644
index 3ce180d719b700a18ce7518e53e113b189265cfc..0000000000000000000000000000000000000000
--- a/src/runner_doiact_rt.c
+++ /dev/null
@@ -1,40 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
- * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- * 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
- * 2020 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-
-/* Config parameters. */
-#include "../config.h"
-
-/* Local headers. */
-#include "active.h"
-#include "cell.h"
-#include "engine.h"
-#include "rt.h"
-#include "runner.h"
-#include "space_getsid.h"
-#include "timers.h"
-
-/* Import the rt injection loop functions. */
-#define FUNCTION inject
-#define FUNCTION_TASK_LOOP TASK_LOOP_RT_INJECT
-#include "runner_doiact_functions_rt.h"
-#undef FUNCTION
-#undef FUNCTION_TASK_LOOP
diff --git a/src/runner_doiact_rt.h b/src/runner_doiact_rt.h
deleted file mode 100644
index 3762b1b2af7cd9a184741027be57e2912eaa5d5b..0000000000000000000000000000000000000000
--- a/src/runner_doiact_rt.h
+++ /dev/null
@@ -1,71 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- * 2018 Loic Hausammann (loic.hausammann@epfl.ch)
- * 2020 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-
-/* Before including this file, define FUNCTION, which is the
- name of the interaction function. This creates the interaction functions
- runner_dopair_FUNCTION, runner_doself_FUNCTION and runner_dosub_FUNCTION
- calling the pairwise interaction function runner_iact_FUNCTION. */
-
-#define PASTE(x, y) x##_##y
-
-#define _DOSELF1_RT(f) PASTE(runner_doself_rt, f)
-#define DOSELF1_RT _DOSELF1_RT(FUNCTION)
-
-#define _DOPAIR1_RT(f) PASTE(runner_dopair_rt, f)
-#define DOPAIR1_RT _DOPAIR1_RT(FUNCTION)
-
-#define _DOPAIR1_NONSYM_RT(f) PASTE(runner_dopair_nonsym_rt, f)
-#define DOPAIR1_NONSYM_RT _DOPAIR1_NONSYM_RT(FUNCTION)
-
-#define _DOSELF1_BRANCH_RT(f) PASTE(runner_doself_branch_rt, f)
-#define DOSELF1_BRANCH_RT _DOSELF1_BRANCH_RT(FUNCTION)
-
-#define _DOPAIR1_BRANCH_RT(f) PASTE(runner_dopair_branch_rt, f)
-#define DOPAIR1_BRANCH_RT _DOPAIR1_BRANCH_RT(FUNCTION)
-
-#define _DOSUB_PAIR1_RT(f) PASTE(runner_dosub_pair_rt, f)
-#define DOSUB_PAIR1_RT _DOSUB_PAIR1_RT(FUNCTION)
-
-#define _DOSUB_SELF1_RT(f) PASTE(runner_dosub_self_rt, f)
-#define DOSUB_SELF1_RT _DOSUB_SELF1_RT(FUNCTION)
-
-#define _TIMER_DOSELF_RT(f) PASTE(timer_doself_rt, f)
-#define TIMER_DOSELF_RT _TIMER_DOSELF_RT(FUNCTION)
-
-#define _TIMER_DOPAIR_RT(f) PASTE(timer_dopair_rt, f)
-#define TIMER_DOPAIR_RT _TIMER_DOPAIR_RT(FUNCTION)
-
-#define _TIMER_DOSUB_SELF_RT(f) PASTE(timer_dosub_self_rt, f)
-#define TIMER_DOSUB_SELF_RT _TIMER_DOSUB_SELF_RT(FUNCTION)
-
-#define _TIMER_DOSUB_PAIR_RT(f) PASTE(timer_dosub_pair_rt, f)
-#define TIMER_DOSUB_PAIR_RT _TIMER_DOSUB_PAIR_RT(FUNCTION)
-
-#define _IACT_RT(f) PASTE(runner_iact_rt, f)
-#define IACT_RT _IACT_RT(FUNCTION)
-
-void DOSELF1_BRANCH_RT(struct runner *r, struct cell *c, int timer);
-void DOPAIR1_BRANCH_RT(struct runner *r, struct cell *ci, struct cell *cj,
- int timer);
-
-void DOSUB_SELF1_RT(struct runner *r, struct cell *ci, int timer);
-void DOSUB_PAIR1_RT(struct runner *r, struct cell *ci, struct cell *cj,
- int timer);
diff --git a/src/runner_ghost.c b/src/runner_ghost.c
index 809170d7c83237437ef55ba50c70427626c986a8..99e1d3df8d48d80b9e71d460ce6a519e58f370c6 100644
--- a/src/runner_ghost.c
+++ b/src/runner_ghost.c
@@ -131,7 +131,8 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
if ((right = (float *)malloc(sizeof(float) * c->stars.count)) == NULL)
error("Can't allocate memory for right.");
for (int k = 0; k < c->stars.count; k++)
- if (spart_is_active(&sparts[k], e) && feedback_is_active(&sparts[k], e)) {
+ if (spart_is_active(&sparts[k], e) &&
+ (feedback_is_active(&sparts[k], e) || with_rt)) {
sid[scount] = k;
h_0[scount] = sparts[k].h;
left[scount] = 0.f;
@@ -156,6 +157,8 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
/* Is this part within the timestep? */
if (!spart_is_active(sp, e))
error("Ghost applied to inactive particle");
+ if (!feedback_is_active(sp, e) && !with_rt)
+ error("Ghost applied to particle inactive for feedback and RT");
#endif
/* Get some useful values */
@@ -252,7 +255,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
feedback_reset_feedback(sp, feedback_props);
}
- if (with_rt && !rt_props->hydro_controlled_injection) {
+ if (with_rt) {
rt_reset_spart(sp);
@@ -423,7 +426,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
feedback_reset_feedback(sp, feedback_props);
}
- if (with_rt && !rt_props->hydro_controlled_injection) {
+ if (with_rt) {
rt_reset_spart(sp);
diff --git a/src/runner_main.c b/src/runner_main.c
index 52f5629c6bd9df8a6307174ede7666ad987c5cf9..2cf729f9a9dc2034a9a0d4eec0242430a72f60cb 100644
--- a/src/runner_main.c
+++ b/src/runner_main.c
@@ -123,13 +123,6 @@
#undef FUNCTION_TASK_LOOP
#undef FUNCTION
-/* Import radiative transfer loop functions. */
-#define FUNCTION inject
-#define FUNCTION_TASK_LOOP TASK_LOOP_RT_INJECT
-#include "runner_doiact_rt.h"
-#undef FUNCTION
-#undef FUNCTION_TASK_LOOP
-
/* Import the RT gradient loop functions */
#define FUNCTION rt_gradient
#define FUNCTION_TASK_LOOP TASK_LOOP_RT_GRADIENT
@@ -269,8 +262,6 @@ void *runner_main(void *data) {
runner_do_bh_swallow_self(r, ci, 1);
else if (t->subtype == task_subtype_bh_feedback)
runner_doself_branch_bh_feedback(r, ci);
- else if (t->subtype == task_subtype_rt_inject)
- runner_doself_branch_rt_inject(r, ci, 1);
else if (t->subtype == task_subtype_rt_gradient)
runner_doself1_branch_rt_gradient(r, ci);
else if (t->subtype == task_subtype_rt_transport)
@@ -319,8 +310,6 @@ void *runner_main(void *data) {
runner_do_bh_swallow_pair(r, ci, cj, 1);
else if (t->subtype == task_subtype_bh_feedback)
runner_dopair_branch_bh_feedback(r, ci, cj);
- else if (t->subtype == task_subtype_rt_inject)
- runner_dopair_branch_rt_inject(r, ci, cj, 1);
else if (t->subtype == task_subtype_rt_gradient)
runner_dopair1_branch_rt_gradient(r, ci, cj);
else if (t->subtype == task_subtype_rt_transport)
@@ -367,8 +356,6 @@ void *runner_main(void *data) {
runner_do_bh_swallow_self(r, ci, 1);
else if (t->subtype == task_subtype_bh_feedback)
runner_dosub_self_bh_feedback(r, ci, 1);
- else if (t->subtype == task_subtype_rt_inject)
- runner_dosub_self_rt_inject(r, ci, 1);
else if (t->subtype == task_subtype_rt_gradient)
runner_dosub_self1_rt_gradient(r, ci, 1);
else if (t->subtype == task_subtype_rt_transport)
@@ -415,8 +402,6 @@ void *runner_main(void *data) {
runner_do_bh_swallow_pair(r, ci, cj, 1);
else if (t->subtype == task_subtype_bh_feedback)
runner_dosub_pair_bh_feedback(r, ci, cj, 1);
- else if (t->subtype == task_subtype_rt_inject)
- runner_dosub_pair_rt_inject(r, ci, cj, 1);
else if (t->subtype == task_subtype_rt_gradient)
runner_dosub_pair1_rt_gradient(r, ci, cj, 1);
else if (t->subtype == task_subtype_rt_transport)
diff --git a/src/scheduler.c b/src/scheduler.c
index bb93cd40cece2541ece1d7102c47d4f5aa2e4437..530949781215d0a0ba9d78d878c40f0cad84116e 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -1558,8 +1558,6 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
t->subtype == task_subtype_force ||
t->subtype == task_subtype_limiter)
cost = 1.f * (wscale * count_i) * count_i;
- else if (t->subtype == task_subtype_rt_inject)
- cost = 1.f * wscale * scount_i * count_i;
else if (t->subtype == task_subtype_rt_gradient)
cost = 1.f * wscale * count_i * count_i;
else if (t->subtype == task_subtype_rt_transport)
@@ -1637,8 +1635,6 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
else
cost = 2.f * (wscale * count_i) * count_j * sid_scale[t->flags];
- } else if (t->subtype == task_subtype_rt_inject) {
- cost = 1.f * wscale * scount_i * count_j;
} else if (t->subtype == task_subtype_rt_gradient) {
cost = 1.f * wscale * count_i * count_j;
} else if (t->subtype == task_subtype_rt_transport) {
@@ -1719,8 +1715,6 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
} else {
cost = 2.f * (wscale * count_i) * count_j * sid_scale[t->flags];
}
- } else if (t->subtype == task_subtype_rt_inject) {
- cost = 1.f * wscale * scount_i * count_j;
} else if (t->subtype == task_subtype_rt_gradient) {
cost = 1.f * wscale * count_i * count_j;
} else if (t->subtype == task_subtype_rt_transport) {
@@ -1755,8 +1749,6 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
t->subtype == task_subtype_force ||
t->subtype == task_subtype_limiter) {
cost = 1.f * (wscale * count_i) * count_i;
- } else if (t->subtype == task_subtype_rt_inject) {
- cost = 1.f * wscale * scount_i * count_i;
} else if (t->subtype == task_subtype_rt_gradient) {
cost = 1.f * wscale * scount_i * count_i;
} else if (t->subtype == task_subtype_rt_transport) {
diff --git a/src/space.c b/src/space.c
index 804f42526e58aa59092d820e9510c2e605704044..7781e96ca80bab0e891cb9cb366f4d571ad840ae 100644
--- a/src/space.c
+++ b/src/space.c
@@ -754,7 +754,7 @@ void space_synchronize_particle_positions(struct space *s) {
void space_convert_rt_star_quantities_after_zeroth_step_mapper(
void *restrict map_data, int scount, void *restrict extra_data) {
-
+#ifdef SWIFT_RT_DEBUG_CHECKS
struct spart *restrict sparts = (struct spart *)map_data;
const struct engine *restrict e = (struct engine *)extra_data;
const int with_cosmology = (e->policy & engine_policy_cosmology);
@@ -787,11 +787,13 @@ void space_convert_rt_star_quantities_after_zeroth_step_mapper(
e->rt_props, e->physical_constants,
e->internal_units);
}
+#endif
}
void space_convert_rt_hydro_quantities_after_zeroth_step_mapper(
void *restrict map_data, int count, void *restrict extra_data) {
+#ifdef SWIFT_RT_DEBUG_CHECKS
struct part *restrict parts = (struct part *)map_data;
const struct engine *restrict e = (struct engine *)extra_data;
const struct rt_props *restrict rt_props = e->rt_props;
@@ -803,39 +805,26 @@ void space_convert_rt_hydro_quantities_after_zeroth_step_mapper(
if (parts[k].time_bin <= num_time_bins)
rt_init_part_after_zeroth_step(p, rt_props);
}
+#endif
}
/**
* @brief Initializes values of radiative transfer data for particles
* that needs to be set before the first actual step is done, but will
* be reset in forthcoming steps when the corresponding particle is
- * active.
- * In hydro controlled injection, in particular we need the stellar
- * emisison rates to be set from the start, not only after the stellar
- * particle has been active. This function requires that the time bins
- * for star particles have been set already and is called after the
- * zeroth time step.
- * In either star controlled injection or hydro controlled injection,
- * for the debug RT scheme some data fields need to be reset after the
- * zeroth step. In particular the interaction count between stars and
- * hydro particles needs to be reset to zero; Otherwise only the active
- * particles/stars will be reset when called in their respective ghosts,
- * while the others remain nonzero, such that the respective sums over
- * all hydro particles and the sum over all star particles won't be
- * equal any longer.
- * TODO MLADEN: Clean this up once you finish with the hydro/star
- * controlled injection and debugging mode.
+ * active. It used to be necessary for a "hydro controlled injection"
+ * idea where the star time steps were necessary to be known, now it's
+ * only used to set up debugging checks correctly.
*
* @param s The #space.
* @param verbose Are we talkative?
*/
void space_convert_rt_quantities_after_zeroth_step(struct space *s,
int verbose) {
-
- const struct rt_props *rt_props = s->e->rt_props;
+#ifdef SWIFT_RT_DEBUG_CHECKS
const ticks tic = getticks();
- if (s->nr_parts > 0 && rt_props->convert_parts_after_zeroth_step)
+ if (s->nr_parts > 0)
/* Particle loop. Reset hydro particle values so we don't inject too much
* radiation into the gas, and other initialisations after zeroth step. */
threadpool_map(&s->e->threadpool,
@@ -843,7 +832,7 @@ void space_convert_rt_quantities_after_zeroth_step(struct space *s,
s->parts, s->nr_parts, sizeof(struct part),
threadpool_auto_chunk_size, /*extra_data=*/s->e);
- if (s->nr_sparts > 0 && rt_props->convert_stars_after_zeroth_step)
+ if (s->nr_sparts > 0)
/* Star particle loop. Hydro controlled injection requires star particles
* to have their emission rates computed and ready for interactions. */
threadpool_map(&s->e->threadpool,
@@ -854,6 +843,7 @@ void space_convert_rt_quantities_after_zeroth_step(struct space *s,
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
clocks_getunit());
+#endif
}
void space_convert_quantities_mapper(void *restrict map_data, int count,
diff --git a/src/space_recycle.c b/src/space_recycle.c
index 9a46a55a14e8e4002d297a9dd736d119beff5d9c..45f4f0a74232ca30e0f9d47220d5d384158cad20 100644
--- a/src/space_recycle.c
+++ b/src/space_recycle.c
@@ -200,7 +200,6 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->stars.ti_end_min = -1;
c->black_holes.ti_end_min = -1;
c->hydro.rt_in = NULL;
- c->hydro.rt_inject = NULL;
c->hydro.rt_ghost1 = NULL;
c->hydro.rt_gradient = NULL;
c->hydro.rt_ghost2 = NULL;
diff --git a/src/task.c b/src/task.c
index 650a7f43a0a21869a7794fb867db4e764364d0bd..a703602b6ed9c4b0b007c70b879379f51da2c33c 100644
--- a/src/task.c
+++ b/src/task.c
@@ -156,7 +156,6 @@ const char *subtaskID_names[task_subtype_count] = {
"do_bh_swallow",
"bh_feedback",
"sink_merger",
- "rt_inject",
"sink_compute_formation",
"sink_accretion",
"rt_gradient",
@@ -293,10 +292,6 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
case task_subtype_sink_compute_formation:
return task_action_all;
- case task_subtype_rt_inject:
- return task_action_all;
- break;
-
case task_subtype_rt_transport:
case task_subtype_rt_gradient:
return task_action_part;
@@ -594,9 +589,6 @@ void task_unlock(struct task *t) {
cell_unlocktree(ci);
} else if (subtype == task_subtype_do_bh_swallow) {
cell_bunlocktree(ci);
- } else if (subtype == task_subtype_rt_inject) {
- cell_unlocktree(ci);
- cell_sunlocktree(ci);
} else if (subtype == task_subtype_limiter) {
#ifdef SWIFT_TASKS_WITHOUT_ATOMICS
cell_unlocktree(ci);
@@ -651,11 +643,6 @@ void task_unlock(struct task *t) {
} else if (subtype == task_subtype_do_bh_swallow) {
cell_bunlocktree(ci);
cell_bunlocktree(cj);
- } else if (subtype == task_subtype_rt_inject) {
- cell_sunlocktree(ci);
- cell_sunlocktree(cj);
- cell_unlocktree(ci);
- cell_unlocktree(cj);
} else if (subtype == task_subtype_limiter) {
#ifdef SWIFT_TASKS_WITHOUT_ATOMICS
cell_unlocktree(ci);
@@ -869,14 +856,6 @@ int task_lock(struct task *t) {
if (ci->hydro.hold) return 0;
if (cell_locktree(ci) != 0) return 0;
#endif
- } else if (subtype == task_subtype_rt_inject) {
- if (ci->stars.hold) return 0;
- if (ci->hydro.hold) return 0;
- if (cell_slocktree(ci) != 0) return 0;
- if (cell_locktree(ci) != 0) {
- cell_sunlocktree(ci);
- return 0;
- }
} else { /* subtype == hydro */
if (ci->hydro.hold) return 0;
if (cell_locktree(ci) != 0) return 0;
@@ -1033,15 +1012,6 @@ int task_lock(struct task *t) {
cell_bunlocktree(ci);
return 0;
}
- } else if (subtype == task_subtype_rt_inject) {
- /* Lock the stars and the gas particles in both cells */
- if (ci->stars.hold || cj->stars.hold) return 0;
- if (ci->hydro.hold || cj->hydro.hold) return 0;
- if (cell_slocktree(ci) != 0) return 0;
- if (cell_slocktree(cj) != 0) {
- cell_sunlocktree(ci);
- return 0;
- }
if (cell_locktree(ci) != 0) {
cell_sunlocktree(ci);
cell_sunlocktree(cj);
@@ -1257,9 +1227,6 @@ void task_get_group_name(int type, int subtype, char *cluster) {
case task_subtype_bh_feedback:
strcpy(cluster, "BHFeedback");
break;
- case task_subtype_rt_inject:
- strcpy(cluster, "RTinject");
- break;
case task_subtype_rt_gradient:
strcpy(cluster, "RTgradient");
break;
@@ -1856,7 +1823,6 @@ enum task_categories task_get_category(const struct task *t) {
case task_subtype_sink_accretion:
return task_category_sink;
- case task_subtype_rt_inject:
case task_subtype_rt_gradient:
case task_subtype_rt_transport:
return task_category_rt;
diff --git a/src/task.h b/src/task.h
index 2725696f71dea423db2c303297015f62767515df..e3b066fc5f92662344d42ebedfc0a004b09881b8 100644
--- a/src/task.h
+++ b/src/task.h
@@ -152,7 +152,6 @@ enum task_subtypes {
task_subtype_do_bh_swallow,
task_subtype_bh_feedback,
task_subtype_sink_merger,
- task_subtype_rt_inject,
task_subtype_sink_compute_formation,
task_subtype_sink_accretion,
task_subtype_rt_gradient,
diff --git a/src/timers.c b/src/timers.c
index 4fc5c5ffd43924cd98830935804190ad9cfe3d8b..202c36fed284af315cf6e9fe092a7df8fa1c8055 100644
--- a/src/timers.c
+++ b/src/timers.c
@@ -125,10 +125,6 @@ const char* timers_names[timer_count] = {
"fof_self",
"fof_pair",
"drift_sink",
- "doself_rt_inject",
- "dopair_rt_inject",
- "dosub_self_rt_inject",
- "dosub_pair_rt_inject",
"rt_ghost1",
"rt_ghost2",
"doself_rt_gradient",
diff --git a/src/timers.h b/src/timers.h
index 4af7947279beca21109bc01745c7afe392c5e6de..4dcd1faccd30ed03a5f3bd6dacb7a9850a7d6810 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -125,10 +125,6 @@ enum {
timer_fof_self,
timer_fof_pair,
timer_drift_sink,
- timer_doself_rt_inject,
- timer_dopair_rt_inject,
- timer_dosub_self_rt_inject,
- timer_dosub_pair_rt_inject,
timer_do_rt_ghost1,
timer_do_rt_ghost2,
timer_doself_rt_gradient,
diff --git a/tools/task_plots/swift_hardcoded_data.py b/tools/task_plots/swift_hardcoded_data.py
index 8beb97e07a95c3f1367a62f2f51df01ac45aac72..dbdd016e1a1843c587980841afc126e5e0b89c3a 100644
--- a/tools/task_plots/swift_hardcoded_data.py
+++ b/tools/task_plots/swift_hardcoded_data.py
@@ -114,7 +114,6 @@ SUBTYPES = [
"do_bh_swallow",
"bh_feedback",
"sink_merger",
- "rt_inject",
"sink_compute_formation",
"sink_accretion",
"rt_gradient",