diff --git a/.gitignore b/.gitignore
index 2f2b42d95d656022649b3b77ef5c07e914cf21af..8feacfd8ae84bb9dccfd8a021ea8a5314715fc68 100644
--- a/.gitignore
+++ b/.gitignore
@@ -222,6 +222,7 @@ theory/Multipoles/mac_potential.pdf
theory/Cosmology/cosmology.pdf
theory/Cooling/eagle_cooling.pdf
theory/Gizmo/gizmo-implementation-details/gizmo-implementation-details.pdf
+theory/RadiativeTransfer/GEARRT/GEARRT.pdf
m4/libtool.m4
m4/ltoptions.m4
diff --git a/configure.ac b/configure.ac
index 72c6fdafe0e2b6ca89050ae00bfbdc53e14f0bad..bc63e5143b691f355efd0b7cf074da8d40432e20 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1308,7 +1308,6 @@ if test "x$with_grackle" != "xno"; then
have_grackle="yes"
- echo $GRACKLE_LIBS
AS_VAR_APPEND([GRACKLE_LIBS], ["$FCLIBS"])
AC_CHECK_LIB(
@@ -1317,8 +1316,24 @@ if test "x$with_grackle" != "xno"; then
[AC_DEFINE([HAVE_GRACKLE],1,[The GRACKLE library appears to be present.])
AC_DEFINE([CONFIG_BFLOAT_8],1,[Use doubles in grackle])
],
- [AC_MSG_ERROR(Cannot find grackle library! You need to have a version >= 3.1.1)],
+ [AC_MSG_ERROR(Cannot find grackle library! Please consult the documentation for specific version required.)],
[$GRACKLE_LIBS])
+
+ AC_CHECK_LIB(
+ [grackle],
+ [set_velocity_units],
+ [ : ], # : means do nothing. Leaving this argument empty triggers AC default behaviour, which breaks stuff down the line.
+ [AC_MSG_ERROR(Wrong grackle library version. Please consult the documentation for specifics.)],
+ [$GRACKLE_LIBS])
+
+ AC_CHECK_LIB(
+ [grackle],
+ [get_grackle_version],
+ [ : ], # : means do nothing
+ [AC_MSG_ERROR(Wrong grackle library version. Please consult the documentation for specifics.)],
+ [$GRACKLE_LIBS])
+
+
fi
AC_SUBST([GRACKLE_LIBS])
AC_SUBST([GRACKLE_INCS])
diff --git a/doc/RTD/source/RadiativeTransfer/GEAR_RT.rst b/doc/RTD/source/RadiativeTransfer/GEAR_RT.rst
index 7fde36f0968b4b6e3b1515b99d3c5078c3affc83..0fd66e35c882e5d1fe601baedc26f0ac4e8137c5 100644
--- a/doc/RTD/source/RadiativeTransfer/GEAR_RT.rst
+++ b/doc/RTD/source/RadiativeTransfer/GEAR_RT.rst
@@ -30,7 +30,7 @@ Compiling for GEAR RT
you to select a hydro Riemann solver, e.g ``--with-riemann-solver=hllc``.
- The thermochemistry requires the `grackle <https://github.com/grackle-project/grackle>`_
- library version 3.2. Grackle is a chemistry and cooling library presented in
+ library version above 3.2.1. [#f4]_ Grackle is a chemistry and cooling library presented in
`B. Smith et al. 2017 <https://ui.adsabs.harvard.edu/abs/2017MNRAS.466.2217S>`_.
Please note that the current implementation is not (yet) as
advanced as the :ref:`GEAR subgrid model grackle cooling <gear_grackle_cooling>`,
@@ -38,6 +38,17 @@ Compiling for GEAR RT
grackle cooling in combination with GEAR RT. You can however follow the Grackle
installation instructions documented there.
+.. warning::
+ (State 2023) Grackle is experiencing current development, and the API is subject
+ to changes in the future. For convenience, a frozen version is hosted as a fork
+ on github here: https://github.com/mladenivkovic/grackle-swift .
+ The version available there will be tried and tested and ensured to work with
+ GEAR-RT.
+
+ Additionally, that repository hosts files necessary to install that specific
+ version of grackle with spack.
+
+
Compulsory Runtime Parameters
@@ -468,3 +479,10 @@ useful:
`I. Iliev et al. 2006 <https://ui.adsabs.harvard.edu/abs/2006MNRAS.369.1625I>`_
paper, but that is very specialized and shouldn't have much use in real
applications.
+
+.. [#f4] Grackle version 3.2.1 still contained a bug related to the use of their
+ "threadsafe functions" that could lead to disastrous outcomes. That was fixed
+ in commit `a59489f`. So versions after 3.2.1 should work as expected.
+ To be safe, we recommend you use the forked grackle repository specifically
+ intended to freeze a stable version for use with SWIFT. You can find that fork
+ on github: https://github.com/mladenivkovic/grackle-swift .
diff --git a/examples/EAGLE_low_z/EAGLE_12/README b/examples/EAGLE_low_z/EAGLE_12/README
index 2e9a053c20ea6859ea3ce0869797f6aabfe25213..bd685ad14400e3d9154ee83b0b0502d0be4a920f 100644
--- a/examples/EAGLE_low_z/EAGLE_12/README
+++ b/examples/EAGLE_low_z/EAGLE_12/README
@@ -26,6 +26,9 @@ o To use GEAR-RT, configure SWIFT with
--with-stars=basic --with-hydro=gizmo-mfv --with-riemann-solver=hllc --with-rt=GEAR_1 --with-rt-riemann-solver=GLF --with-feedback=none
[technically, any other feedback scheme should work as well.]
+[the .yml file is set up to handle up to 3 photon groups, so
+ you could also be compiling with --with-rt=GEAR_3 if you want
+ to test more photon groups as well.]
o To use the DEBUG RT scheme, configure SWIFT with
diff --git a/examples/EAGLE_low_z/EAGLE_12/eagle_12_rt_test.yml b/examples/EAGLE_low_z/EAGLE_12/eagle_12_rt_test.yml
index bd5c4139fc2a52cd57d41c6c61f7906029262ba6..6fb36b2a7dd566aeb85e2f2976bef75a5dbe81db 100644
--- a/examples/EAGLE_low_z/EAGLE_12/eagle_12_rt_test.yml
+++ b/examples/EAGLE_low_z/EAGLE_12/eagle_12_rt_test.yml
@@ -98,5 +98,6 @@ GEARRT:
stellar_spectrum_type: 1 # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
stellar_spectrum_blackbody_temperature_K: 1.e5 # (Conditional) if stellar_spectrum_type=1, use this temperature (in K) for the blackbody spectrum.
stars_max_timestep: 5.468750e-05 # (Optional) restrict the maximal timestep of stars to this value (in internal units)
+ max_tchem_recursion: 10 # (Optional) if > 0, sets the maximal recursion depth when re-computing the thermochemistry if |u_new/u_old - 1| > 0.1.
diff --git a/examples/EAGLE_low_z/EAGLE_25/README b/examples/EAGLE_low_z/EAGLE_25/README
index 80907ce9ffff04382afceb7c906dca089a05b8bb..5ec0be761afa4255e4ddee535a818291265661a2 100644
--- a/examples/EAGLE_low_z/EAGLE_25/README
+++ b/examples/EAGLE_low_z/EAGLE_25/README
@@ -14,3 +14,24 @@ running these ICs on 64 cores.
MD5 checksum of the ICs:
e6a5de0e962d8ffb7589671b9613daa0 EAGLE_ICs_25.hdf5
+
+
+A second example is set up to run with `./run_rt_test.sh`. It makes use
+of the `eagle_25_rt_test.yml` file, and is intended to stress-test the
+RT implementation rather than produce actual physically meaningful results.
+
+
+o To use GEAR-RT, configure SWIFT with
+
+ --with-stars=basic --with-hydro=gizmo-mfv --with-riemann-solver=hllc --with-rt=GEAR_1 --with-rt-riemann-solver=GLF --with-feedback=none
+
+[technically, any other feedback scheme should work as well.]
+[the .yml file is set up to handle up to 3 photon groups, so
+ you could also be compiling with --with-rt=GEAR_3 if you want
+ to test more photon groups as well.]
+
+
+o To use the DEBUG RT scheme, configure SWIFT with
+ --with-stars=basic --with-rt=debug --with-feedback=none
+
+[technically, any other feedback scheme should work as well.]
diff --git a/examples/EAGLE_low_z/EAGLE_25/eagle_25_rt_test.yml b/examples/EAGLE_low_z/EAGLE_25/eagle_25_rt_test.yml
new file mode 100644
index 0000000000000000000000000000000000000000..ae1cbe4f4f8847fbf07fbe1f14b25ebd3a1373e3
--- /dev/null
+++ b/examples/EAGLE_low_z/EAGLE_25/eagle_25_rt_test.yml
@@ -0,0 +1,103 @@
+MetaData:
+ run_name: EAGLE-L0025N0376-Ref-RT-TEST
+
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.98841e43 # 10^10 M_sun in grams
+ UnitLength_in_cgs: 3.08567758e24 # Mpc in centimeters
+ UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Cosmological parameters
+Cosmology:
+ h: 0.6777 # Reduced Hubble constant
+ a_begin: 0.9090909 # Initial scale-factor of the simulation
+ a_end: 1.0 # Final scale factor of the simulation
+ Omega_cdm: 0.2587481 # Cold Dark Matter density parameter
+ Omega_lambda: 0.693 # Dark-energy density parameter
+ Omega_b: 0.0482519 # Baryon density parameter
+
+# Parameters governing the time integration
+TimeIntegration:
+ max_nr_rt_subcycles: 32
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 1e-2 # The end time of the simulation (in internal units).
+ dt_min: 1e-12 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: eagle # Common part of the name of output files
+ scale_factor_first: 0.91 # Scale-factor of the first snaphot (cosmological run)
+ time_first: 0.01 # Time of the first output (non-cosmological run) (in internal units)
+ delta_time: 1.01 # Time difference between consecutive outputs (in internal units)
+
+Scheduler:
+ links_per_tasks: 500
+ dependency_graph_frequency: 0
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ scale_factor_first: 0.92 # Scale-factor of the first stat dump (cosmological run)
+ time_first: 0.01 # Time of the first stat dump (non-cosmological run) (in internal units)
+ delta_time: 1.05 # Time between statistics output
+
+# Parameters for the self-gravity scheme
+Gravity:
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ MAC: adaptive
+ epsilon_fmm: 0.001
+ theta_cr: 0.7 # Opening angle (Multipole acceptance criterion)
+ mesh_side_length: 256
+ comoving_DM_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_DM_softening: 0.0007 # Max physical DM softening length (in internal units).
+ comoving_baryon_softening: 0.0026994 # Comoving DM softening length (in internal units).
+ max_physical_baryon_softening: 0.0007 # Max physical DM softening length (in internal units).
+
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ h_min_ratio: 0.1 # Minimal smoothing length in units of softening.
+ h_max: 0.5 # Maximal smoothing length in co-moving internal units.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 100 # (internal units)
+ particle_splitting: 0
+ particle_splitting_mass_threshold: 7e-4 # Internal units (i.e. 7e6 Msun ~ 4 times the initial gas mass)
+
+# Parameters of the stars neighbour search
+Stars:
+ resolution_eta: 1.1642 # Target smoothing length in units of the mean inter-particle separation
+ h_tolerance: 7e-3
+ overwrite_birth_time: 1
+ birth_time: 0.33333 # Pretend all the stars were born at z = 2
+
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./EAGLE_ICs_25.hdf5 # The file to read
+ periodic: 1
+ cleanup_h_factors: 1 # Remove the h-factors inherited from Gadget
+ cleanup_velocity_factors: 1 # Remove the sqrt(a) factor in the velocities inherited from Gadget
+
+GEARRT:
+ f_reduce_c: 1e-2 # reduce the speed of light for the RT solver by multiplying c with this factor
+ CFL_condition: 0.4 # CFL condition for RT, independent of hydro
+ f_limit_cooling_time: 0.01 # (Optional) multiply the cooling time by this factor when estimating maximal next time step. Set to 0.0 to turn computation of cooling time off.
+ photon_groups_Hz: [3.288e15, 5.944e15, 13.157e15] # Lower photon frequency group bin edges in Hz. Needs to have exactly N elements, where N is the configured number of bins --with-RT=GEAR_N
+ stellar_luminosity_model: const # Which model to use to determine the stellar luminosities.
+ const_stellar_luminosities_LSol: [1.764e+04, 3.631e+04, 8.037e+03 ] # (Conditional) constant star emission rates for each photon frequency group to use if use_constant_emission_rates is set, in units of Solar Luminosity.
+ hydrogen_mass_fraction: 0.75 # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
+ set_equilibrium_initial_ionization_mass_fractions: 0 # (Optional) set the initial ionization fractions depending on gas temperature assuming ionization equilibrium.
+ set_initial_ionization_mass_fractions: 1 # (Optional) manually overwrite initial mass fraction of each species (using the values you set below)
+ mass_fraction_HI: 0.75 # (Conditional) If set_initial_ionization_fractions=1, needed to set initial HI mass fractions to this value
+ mass_fraction_HII: 1.e-6 # (Conditional) If set_initial_ionization_fractions=1, needed to set initial HII mass fractions to this value
+ mass_fraction_HeI: 0.25 # (Conditional) If set_initial_ionization_fractions=1, needed to set initial HeI mass fractions to this value
+ mass_fraction_HeII: 1.e-6 # (Conditional) If set_initial_ionization_fractions=1, needed to set initial HeII mass fractions to this value
+ mass_fraction_HeIII: 1.e-6 # (Conditional) If set_initial_ionization_fractions=1, needed to set initial HeIII mass fractions to this value
+ stellar_spectrum_type: 1 # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
+ stellar_spectrum_blackbody_temperature_K: 1.e5 # (Conditional) if stellar_spectrum_type=1, use this temperature (in K) for the blackbody spectrum.
+ stars_max_timestep: -1. # (Optional) restrict the maximal timestep of stars to this value (in internal units). Set to negative to turn off.
+ max_tchem_recursion: 10 # (Optional) if > 0, sets the maximal recursion depth when re-computing the thermochemistry if |u_new/u_old - 1| > 0.1.
+
diff --git a/examples/EAGLE_low_z/EAGLE_25/run_rt_test.sh b/examples/EAGLE_low_z/EAGLE_25/run_rt_test.sh
new file mode 100755
index 0000000000000000000000000000000000000000..cd5417c4b636747e1e6555ef74c05b4b16c3409f
--- /dev/null
+++ b/examples/EAGLE_low_z/EAGLE_25/run_rt_test.sh
@@ -0,0 +1,27 @@
+#!/bin/bash
+
+# Run a RT test with EAGLE ICs.
+#
+# To use GEAR-RT, configure SWIFT with
+#
+# --with-stars=basic --with-hydro=gizmo-mfv --with-riemann-solver=hllc --with-rt=GEAR_1 --with-rt-riemann-solver=GLF --with-feedback=none
+# [technically, any other feedback scheme should work as well.]
+#
+#
+# To use the DEBUG RT scheme, configure SWIFT with
+# --with-stars=basic --with-rt=debug --with-feedback=none
+# [technically, any other feedback scheme should work as well.]
+
+
+# Generate the initial conditions if they are not present.
+if [ ! -e EAGLE_ICs_25.hdf5 ]
+then
+ echo "Fetching initial conditions for the EAGLE 12Mpc example..."
+ ./getIC.sh
+fi
+
+../../../swift \
+ --hydro --threads=16 --stars --self-gravity \
+ --feedback --radiation \
+ eagle_25_rt_test.yml
+
diff --git a/examples/RadiativeTransferTests/CoolingTest/plotSolution.py b/examples/RadiativeTransferTests/CoolingTest/plotSolution.py
index b8622ff805835589115dac07804ea31095709c7e..5ea851990a6d0284839b723ccb4467da1e52a8db 100755
--- a/examples/RadiativeTransferTests/CoolingTest/plotSolution.py
+++ b/examples/RadiativeTransferTests/CoolingTest/plotSolution.py
@@ -48,7 +48,7 @@ mass_units = unyt.Msun
def mean_molecular_weight(XH0, XHp, XHe0, XHep, XHepp):
"""
- Determines the mean molecular weight for given
+ Determines the mean molecular weight for given
mass fractions of
hydrogen: XH0
H+: XHp
@@ -75,7 +75,7 @@ def mean_molecular_weight(XH0, XHp, XHe0, XHep, XHepp):
def gas_temperature(u, mu, gamma):
"""
- Compute the gas temperature given the specific internal
+ Compute the gas temperature given the specific internal
energy u and the mean molecular weight mu
"""
@@ -89,7 +89,7 @@ def gas_temperature(u, mu, gamma):
def get_snapshot_list(snapshot_basename="output"):
"""
- Find the snapshot(s) that are to be plotted
+ Find the snapshot(s) that are to be plotted
and return their names as list
"""
@@ -175,7 +175,7 @@ def get_snapshot_data(snaplist):
Returns:
numpy arrays of:
time
- temperatures
+ temperatures
mean molecular weights
mass fractions
"""
@@ -189,7 +189,7 @@ def get_snapshot_data(snaplist):
# allow to read in solutions with only cooling, without RT
with_rt = False
- times = np.zeros(nsnaps) * mass_units
+ times = np.zeros(nsnaps) * unyt.Myr
temperatures = np.zeros(nsnaps) * unyt.K
mean_molecular_weights = np.zeros(nsnaps)
mass_fractions = np.zeros((nsnaps, 5))
diff --git a/examples/RadiativeTransferTests/RandomizedBox_3D/README b/examples/RadiativeTransferTests/RandomizedBox_3D/README
index 1d8fb2855841811b5e7d05ad521e4f5af5e0a485..28fa043ca315183056c2f454f2e29c21f0cdaaf3 100644
--- a/examples/RadiativeTransferTests/RandomizedBox_3D/README
+++ b/examples/RadiativeTransferTests/RandomizedBox_3D/README
@@ -23,4 +23,6 @@ Some specific notes for debugging:
- when running on 2 MPI ranks (at least in the first few steps) cells 27 and
151 have proxies on the other rank. However, those are made for the gravity
tasks, not for the RT tasks. So their foreign counterparts will have no RT
- tasks at all.
+ tasks at all. (run with debug-RT scheme and 4 sub-cycles. To catch the error,
+ remove additional check in cell_is_rt_active() to see whether cell actually
+ does RT.)
diff --git a/m4/ax_gcc_archflag.m4 b/m4/ax_gcc_archflag.m4
index 2754d8a40d912e2674719cea5e1338c6ab0dd693..342396d6bbd2f79bdc916315a5233e54e8d5d44f 100644
--- a/m4/ax_gcc_archflag.m4
+++ b/m4/ax_gcc_archflag.m4
@@ -147,6 +147,8 @@ case $host_cpu in
*06??f??:*:*:*|6??f??:*:*:*) ax_gcc_arch="bdver3 bdver2 bdver1 amdfam10 k8" ;;
*070?f??:*:*:*|70?f??:*:*:*) ax_gcc_arch="btver2 btver1 amdfam10 k8" ;;
83?f??:*:*:*) ax_gcc_arch="znver2 znver1 btver2 btver1 amdfam10 k8" ;;
+ a0?f??:*:*:*) ax_gcc_arch="znver3 znver2 znver1 btver2 btver1 amdfam10 k8" ;;
+ a1?f??:*:*:*) ax_gcc_arch="znver4 znver3 znver2 znver1 btver2 btver1 amdfam10 k8" ;;
*???f??:*:*:*) ax_gcc_arch="amdfam10 k8" ;;
esac ;;
*:746e6543:736c7561:48727561) # IDT / VIA (Centaur)
diff --git a/src/active.h b/src/active.h
index 3150a13a4e426ac6240e1547aeedbecdc5a2a874..ae2c2423907f2e94b2e5073307efff5e228f21f9 100644
--- a/src/active.h
+++ b/src/active.h
@@ -199,18 +199,7 @@ __attribute__((always_inline)) INLINE static int cell_is_rt_active(
struct cell *c, const struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
- /* Each cell doing RT needs to have the rt_advance_cell_time task.
- * If it doesn't, it's not doing RT. This can happen for foreign
- * cells which have been sent to a foreign node for other interactions,
- * e.g. gravity. However, foreign cells may have tasks on levels below
- * the rt_advance_cell_time, so allow for that exception in this check. */
-
- int has_rt_advance_cell_time = 0;
- if (c->super != NULL)
- has_rt_advance_cell_time = c->super->rt.rt_advance_cell_time != NULL;
-
- if (has_rt_advance_cell_time &&
- (c->rt.ti_rt_end_min < e->ti_current_subcycle)) {
+ if (c->rt.ti_rt_end_min < e->ti_current_subcycle) {
error(
"cell %lld in an impossible time-zone! c->ti_rt_end_min=%lld (t=%e) "
"and e->ti_current=%lld (t=%e, a=%e) c->nodeID=%d ACT=%d count=%d",
diff --git a/src/cell_unskip.c b/src/cell_unskip.c
index c4a52110a0caf6927e86f2efe771879462ca124b..406505fd8c4ffd66c7f19b0f3331d0f77b29e5e5 100644
--- a/src/cell_unskip.c
+++ b/src/cell_unskip.c
@@ -3321,12 +3321,19 @@ int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s,
}
/* The rt_advance_cell_time tasks also run on foreign cells */
- if (c->super != NULL && c->super->rt.rt_advance_cell_time != NULL)
+ if (c->super != NULL && c->super->rt.rt_advance_cell_time != NULL) {
scheduler_activate(s, c->super->rt.rt_advance_cell_time);
- /* The rt_collect_times tasks replace the timestep_collect tasks
- * during sub-cycles, so we only activate it when sub-cycling. */
- if (c->rt.rt_collect_times != NULL && sub_cycle)
- scheduler_activate(s, c->rt.rt_collect_times);
+ }
+ if (sub_cycle) {
+ /* The rt_collect_times tasks replace the timestep_collect tasks
+ * during sub-cycles, so we only activate it when sub-cycling. */
+ if (c->top->rt.rt_collect_times != NULL)
+ scheduler_activate(s, c->top->rt.rt_collect_times);
+ } else {
+ /* Otherwise, make sure timestep_collect and timestep is active. */
+ if (c->top->timestep_collect != NULL)
+ scheduler_activate(s, c->top->timestep_collect);
+ }
}
return rebuild;
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index 68628c7881e65b721af8cc6cdad6b297e0bedffb..e205cf1f1a74371d7b117868095328b071bbb51f 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -626,38 +626,17 @@ void engine_addtasks_recv_hydro(
s, task_type_recv, task_subtype_rt_transport, c->mpi.tag, 0, c, NULL);
/* Also create the rt_advance_cell_time tasks for the foreign cells
* for the sub-cycling. */
+
#ifdef SWIFT_RT_DEBUG_CHECKS
if (c->super == NULL)
error("trying to add rt_advance_cell_time above super level...");
-#endif
-
- /* Create the RT collect times task at the top level, if it hasn't
- * already. */
if (c->top->rt.rt_collect_times == NULL) {
- c->top->rt.rt_collect_times =
- scheduler_addtask(s, task_type_rt_collect_times, task_subtype_none,
- 0, 0, c->top, NULL);
+ error("rt_collect_times should exist already");
}
-
- /* Create the RT advance times task at the super level, if it hasn't
- * already. Also set all the dependencies */
if (c->super->rt.rt_advance_cell_time == NULL) {
-
- c->super->rt.rt_advance_cell_time =
- scheduler_addtask(s, task_type_rt_advance_cell_time,
- task_subtype_none, 0, 0, c->super, NULL);
-
- /* Don't run collect times before you run advance cell time */
- scheduler_addunlock(s, c->super->rt.rt_advance_cell_time,
- c->top->rt.rt_collect_times);
-
- /* In normal steps, tend mustn't run before rt_advance_cell_time or the
- * cell's ti_rt_end_min will be updated wrongly. In sub-cycles, we don't
- * have the tend tasks, so there's no worry about that. (Them missing is
- * the reason we need the rt_advanced_cell_time to complete the
- * sub-cycles in the first place) */
- scheduler_addunlock(s, c->super->rt.rt_advance_cell_time, tend);
+ error("rt_advance_cell_times should exist already");
}
+#endif
/* Make sure we sort after receiving RT data. The hydro sorts may or may
* not be active. Blocking them with dependencies deadlocks with MPI. So
@@ -828,6 +807,72 @@ void engine_addtasks_recv_hydro(
#endif
}
+/**
+ * @brief Add time rt_advance_cell_time tasks to super levels of
+ * foreign cells. This function recurses down to the super level
+ * and creates the required tasks, and adds a dependency between
+ * rt_advance_cell_time, rt_collect_times, and tend tasks.
+ *
+ * In normal steps, tend mustn't run before rt_advance_cell_time or the
+ * cell's ti_rt_end_min will be updated wrongly. In sub-cycles, we don't
+ * have the tend tasks, so there's no worry about that. (Them missing is
+ * the reason we need the rt_advanced_cell_time to complete the
+ * sub-cycles in the first place)
+ *
+ * @param e The #engine.
+ * @param c The foreign #cell.
+ * @param tend The top-level time-step communication #task.
+ */
+
+void engine_addtasks_recv_rt_advance_cell_time(struct engine *e, struct cell *c,
+ struct task *const tend) {
+
+#ifdef WITH_MPI
+ struct scheduler *s = &e->sched;
+
+ /* Early abort (are we below the level where tasks are)? */
+ if (!cell_get_flag(c, cell_flag_has_tasks)) return;
+
+ /* Have we reached the super level? */
+ if (c->super == c) {
+
+#ifdef SWIFT_RT_DEBUG_CHECKS
+ if (c->super == NULL)
+ error("trying to add rt_advance_cell_time above super level...");
+ if (c->top->rt.rt_collect_times == NULL)
+ error("rt_collect_times should have been created already????");
+#endif
+
+ /* Create the rt advance times task at the super level, if it hasn't
+ * already. also set all the dependencies */
+ if (c->rt.rt_advance_cell_time == NULL) {
+
+ c->rt.rt_advance_cell_time = scheduler_addtask(
+ s, task_type_rt_advance_cell_time, task_subtype_none, 0, 0, c, NULL);
+
+ /* don't run collect times before you run advance cell time */
+ scheduler_addunlock(s, c->rt.rt_advance_cell_time,
+ c->top->rt.rt_collect_times);
+
+ /* Add the dependency */
+ scheduler_addunlock(s, c->super->rt.rt_advance_cell_time, tend);
+ }
+
+ /* we're done. */
+ return;
+ }
+
+ /* Recurse? */
+ if (c->split)
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL)
+ engine_addtasks_recv_rt_advance_cell_time(e, c->progeny[k], tend);
+
+#else
+ error("SWIFT was not compiled with MPI support.");
+#endif
+}
+
/**
* @brief Add recv tasks for stars pairs to a hierarchy of cells.
*
@@ -1684,7 +1729,7 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
/* In cases where nothing but RT is active, don't allow the timestep
* collect to run before we've finished */
- scheduler_addunlock(s, c->rt.rt_out, c->super->timestep_collect);
+ scheduler_addunlock(s, c->rt.rt_out, c->top->timestep_collect);
/* non-implicit ghost 1 */
c->rt.rt_ghost1 = scheduler_addtask(s, task_type_rt_ghost1,
@@ -4186,6 +4231,57 @@ struct cell_type_pair {
int type;
};
+/**
+ * @brief Recurse down to the super level and add a dependency between
+ * rt_advance_cell_time and tend tasks. Note: This function is intended
+ * for the sending side, i.e. for local cells.
+ *
+ * If we're running with RT subcycling, we need to ensure that nothing
+ * is sent before the advance cell time task has finished. This may
+ * overwrite the correct cell times, particularly so when we're sending
+ * over data for non-RT tasks, e.g. for gravity pair tasks. Therefore the
+ * send/tend task needs to be unlocked by the rt_advance_cell_time task.
+ *
+ * The send/tend task is on the top level, while the rt_advance_cell_time
+ * task is on the super level. This function simply recurses down to the
+ * super level and adds the required dependency.
+ *
+ * @param c cell to check/recurse into
+ * @param tend the send/tend task that needs to be unlocked.
+ * @param e the engine
+ */
+void engine_addunlock_rt_advance_cell_time_tend(struct cell *c,
+ struct task *tend,
+ struct engine *e) {
+
+ /* safety measure */
+ if (!cell_get_flag(c, cell_flag_has_tasks)) return;
+ if (cell_is_empty(c)) return;
+
+ if (c->super == c) {
+ /* Found the super level cell. Add dependency from rt_advance_cell_time, if
+ * it exists. */
+ if (c->super->rt.rt_advance_cell_time != NULL) {
+ scheduler_addunlock(&e->sched, c->super->rt.rt_advance_cell_time, tend);
+ }
+#ifdef SWIFT_RT_DEBUG_CHECKS
+ else {
+ error("Got local super cell without rt_advance_cell_time task");
+ }
+#endif
+
+ } else {
+ /* descend the tree until you find the super level */
+ if (c->split) {
+ for (int k = 0; k < 8; k++) {
+ if (c->progeny[k] != NULL) {
+ engine_addunlock_rt_advance_cell_time_tend(c->progeny[k], tend, e);
+ }
+ }
+ }
+ }
+}
+
void engine_addtasks_send_mapper(void *map_data, int num_elements,
void *extra_data) {
@@ -4217,21 +4313,8 @@ void engine_addtasks_send_mapper(void *map_data, int num_elements,
scheduler_addunlock(&e->sched, ci->timestep_collect, tend);
engine_addlink(e, &ci->mpi.send, tend);
- if (with_rt && (type & proxy_cell_type_hydro)) {
-
- /* If we're running with RT subcycling, we need to ensure that nothing
- * is sent before the advance cell time task has finished. This may
- * overwrite the correct cell times, particularly so when we're sending
- * over data for non-RT tasks, e.g. for gravity pair tasks. */
- if (ci->super->rt.rt_advance_cell_time != NULL) {
- scheduler_addunlock(&e->sched, ci->super->rt.rt_advance_cell_time,
- tend);
-#ifdef SWIFT_RT_DEBUG_CHECKS
- } else {
- error("Got local super cell without rt_advance_cell_time task");
-#endif
- }
- }
+ if (with_rt && (type & proxy_cell_type_hydro))
+ engine_addunlock_rt_advance_cell_time_tend(ci, tend, e);
}
#endif
@@ -4299,6 +4382,35 @@ void engine_addtasks_recv_mapper(void *map_data, int num_elements,
tend = scheduler_addtask(&e->sched, task_type_recv, task_subtype_tend,
ci->mpi.tag, 0, ci, NULL);
engine_addlink(e, &ci->mpi.recv, tend);
+
+ /* If we're running with RT, there may be foreign cells that
+ * don't receive any actual hydro particles. These cells however
+ * still need to have the "advance_cell_time" and "rt_collect_times"
+ * tasks in order for their time variables to be correct, especially
+ * during syb-cycles, where the cell times aren't communicated at the
+ * end of the step. So we create them now. */
+ if (with_rt) {
+#ifdef SWIFT_RT_DEBUG_CHECKS
+ if (ci->top == NULL) error("Working on a cell with top == NULL??");
+#endif
+
+ /* Create the RT collect times task at the top level, if it hasn't
+ * already. */
+ if (ci->top->rt.rt_collect_times == NULL) {
+ ci->top->rt.rt_collect_times =
+ scheduler_addtask(&e->sched, task_type_rt_collect_times,
+ task_subtype_none, 0, 0, ci->top, NULL);
+ }
+ /* We don't need rt_collect_times -> tend dependencies. They never
+ * run at the same time. rt_collect_times runs in sub-cycles,
+ * tend runs on normal steps. */
+
+ /* Make sure the timestep task replacements, i.e. rt_advance_cell_time,
+ * exists on the super levels regardless of proxy type.
+ * This needs to be done before engine_addtasks_recv_hydro so we
+ * can set appropriate unlocks there without re-creating tasks. */
+ engine_addtasks_recv_rt_advance_cell_time(e, ci, tend);
+ }
}
#endif
diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c
index 1919b1b501cf2138d1e687edb98bf63060a8e563..a5e6ad5b2ef6697938340b5b4c3c833d28876694 100644
--- a/src/engine_marktasks.c
+++ b/src/engine_marktasks.c
@@ -1564,8 +1564,8 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
if (cell_is_active_black_holes(t->ci, e)) scheduler_activate(s, t);
}
- /* Time-step or time-step collection? */
- else if (t_type == task_type_timestep || t_type == task_type_collect) {
+ /* Time-step collection? */
+ else if (t_type == task_type_timestep) {
t->ci->hydro.updated = 0;
t->ci->grav.updated = 0;
t->ci->stars.updated = 0;
@@ -1581,6 +1581,25 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
}
+ /* Time-step collection? */
+ else if (t_type == task_type_collect) {
+ t->ci->hydro.updated = 0;
+ t->ci->grav.updated = 0;
+ t->ci->stars.updated = 0;
+ t->ci->sinks.updated = 0;
+ t->ci->black_holes.updated = 0;
+ t->ci->rt.updated = 0; /* this is different from time-step */
+
+ if (!cell_is_empty(t->ci)) {
+ if (cell_is_active_hydro(t->ci, e) ||
+ cell_is_active_gravity(t->ci, e) ||
+ cell_is_active_stars(t->ci, e) || cell_is_active_sinks(t->ci, e) ||
+ cell_is_active_black_holes(t->ci, e) || cell_is_rt_active(t->ci, e))
+ /* this is different from time-step ----^*/
+ scheduler_activate(s, t);
+ }
+ }
+
else if ((t_type == task_type_send && t_subtype == task_subtype_tend) ||
(t_type == task_type_recv && t_subtype == task_subtype_tend)) {
if (!cell_is_empty(t->ci)) {
diff --git a/src/equation_of_state/planetary/hm80.h b/src/equation_of_state/planetary/hm80.h
index 1aa746a76c3c9a03596cceb822e7954621a8a76a..d969e86e88d5807e514fad8606d4a63188e99a5d 100644
--- a/src/equation_of_state/planetary/hm80.h
+++ b/src/equation_of_state/planetary/hm80.h
@@ -42,7 +42,7 @@
// Hubbard & MacFarlane (1980) parameters
struct HM80_params {
float *table_log_P_rho_u;
- int date, num_rho, num_u;
+ int version_date, num_rho, num_u;
float log_rho_min, log_rho_max, log_rho_step, inv_log_rho_step, log_u_min,
log_u_max, log_u_step, inv_log_u_step, bulk_mod, P_min_for_c_min;
enum eos_planetary_material_id mat_id;
@@ -54,29 +54,29 @@ INLINE static void set_HM80_HHe(struct HM80_params *mat,
mat->mat_id = mat_id;
mat->bulk_mod = 0.f;
mat->P_min_for_c_min = 1e3f;
- mat->date = 20201003;
+ mat->version_date = 20230710;
}
INLINE static void set_HM80_ice(struct HM80_params *mat,
enum eos_planetary_material_id mat_id) {
mat->mat_id = mat_id;
mat->bulk_mod = 2.0e9f;
mat->P_min_for_c_min = 0.f;
- mat->date = 20201003;
+ mat->version_date = 20230710;
}
INLINE static void set_HM80_rock(struct HM80_params *mat,
enum eos_planetary_material_id mat_id) {
mat->mat_id = mat_id;
mat->bulk_mod = 3.49e10f;
mat->P_min_for_c_min = 0.f;
- mat->date = 20201003;
+ mat->version_date = 20230710;
}
// Read the table from file
INLINE static void load_table_HM80(struct HM80_params *mat, char *table_file) {
/* File contents:
- header (five lines)
- date
+ header (11 lines)
+ version_date
log_rho_min log_rho_max num_rho log_u_min log_u_max num_u (SI)
P_0_0 P_0_1 ... P_0_num_u # Array of pressures (Pa)
P_1_0 ... ... P_1_num_u
@@ -94,21 +94,21 @@ INLINE static void load_table_HM80(struct HM80_params *mat, char *table_file) {
// Ignore header lines
char buffer[100];
- for (int i = 0; i < 5; i++) {
+ for (int i = 0; i < 11; i++) {
if (fgets(buffer, 100, f) == NULL)
error("Failed to read the HM80 EoS file header %s", table_file);
}
// Table properties
- int date;
- int c = fscanf(f, "%d", &date);
+ int version_date;
+ int c = fscanf(f, "%d", &version_date);
if (c != 1) error("Failed to read the HM80 EoS table %s", table_file);
- if (date != mat->date)
+ if (version_date != mat->version_date)
error(
- "EoS file %s date %d does not match expected %d"
+ "EoS file %s version_date %d does not match expected %d"
"\nPlease download the file using "
"examples/Planetary/EoSTables/get_eos_tables.sh",
- table_file, date, mat->date);
+ table_file, version_date, mat->version_date);
c = fscanf(f, "%f %f %d %f %f %d", &mat->log_rho_min, &mat->log_rho_max,
&mat->num_rho, &mat->log_u_min, &mat->log_u_max, &mat->num_u);
if (c != 6) error("Failed to read the HM80 EoS table %s", table_file);
diff --git a/src/riemann/riemann_exact.h b/src/riemann/riemann_exact.h
index e19ead49941f9e0ca712581af1716f1bd288c0c9..57fe127798a95625bde7887f024c48355f077941 100644
--- a/src/riemann/riemann_exact.h
+++ b/src/riemann/riemann_exact.h
@@ -187,8 +187,8 @@ __attribute__((always_inline)) INLINE static float riemann_guess_p(
*
* @param lower_limit Lower limit for the method (riemann_f(lower_limit) < 0)
* @param upper_limit Upper limit for the method (riemann_f(upper_limit) > 0)
- * @param lowf ??? Bert?
- * @param upf ??? Bert?
+ * @param lowf Function value riemann_f(lower_limit)
+ * @param upf Function value riemann_f(upper_limit)
* @param error_tol Tolerance used to decide if the solution is converged
* @param WL Left state vector
* @param WR Right state vector
@@ -202,15 +202,16 @@ __attribute__((always_inline)) INLINE static float riemann_solve_brent(
float error_tol, const float* WL, const float* WR, float vL, float vR,
float aL, float aR) {
- float a, b, c, d, s;
+ float a, b, c, d, e, s;
float fa, fb, fc, fs;
float tmp, tmp2;
int mflag;
a = lower_limit;
b = upper_limit;
- c = 0.0f;
- d = FLT_MAX;
+ c = 0.0f; // previous value of b: b_{n-1}
+ d = FLT_MAX; // value of b from two iterations ago: b_{n-2}
+ e = FLT_MAX; // previous value of a: a_{n-1}
fa = lowf;
fb = upf;
@@ -243,7 +244,15 @@ __attribute__((always_inline)) INLINE static float riemann_solve_brent(
fc = fa;
mflag = 1;
- while (!(fb == 0.0f) && (fabs(a - b) > error_tol * 0.5f * (a + b))) {
+ /* Loop until convergence, i.e. until an exact zero point is found, or the
+ * interval is sufficiently small, or the interval is unchanged since the
+ * previous iteration */
+ int counter = 0;
+ while ((fb != 0.0f) && (fabs(a - b) > error_tol * 0.5f * (a + b)) &&
+ (a != e || b != c)) {
+ counter++;
+ if (counter > 1000) error("Brent's method did not converge!\n");
+
if ((fa != fc) && (fb != fc)) /* Inverse quadratic interpolation */
s = a * fb * fc / (fa - fb) / (fa - fc) +
b * fa * fc / (fb - fa) / (fb - fc) +
@@ -266,6 +275,7 @@ __attribute__((always_inline)) INLINE static float riemann_solve_brent(
fs = riemann_f(s, WL, WR, vL, vR, aL, aR);
d = c;
c = b;
+ e = a;
fc = fb;
if (fa * fs < 0.) {
b = s;
diff --git a/src/rt/GEAR/rt.h b/src/rt/GEAR/rt.h
index 9d0b21ce119ec36317028f78a9d358dc9cc0207a..72a9e25b221c162063dbffbad0e8d84abe6ed37e 100644
--- a/src/rt/GEAR/rt.h
+++ b/src/rt/GEAR/rt.h
@@ -568,7 +568,7 @@ __attribute__((always_inline)) INLINE static void rt_kick_extra(
/* Update the mass fraction changes due to interparticle fluxes */
const float current_mass = p->conserved.mass;
- if (current_mass <= 0.f || p->rho <= 0.f) {
+ if ((current_mass <= 0.f) || (p->rho <= 0.f)) {
/* Deal with vacuum. Let hydro deal with actuall mass < 0, just do your
* mass fractions thing. */
p->rt_data.tchem.mass_fraction_HI = 0.f;
@@ -683,7 +683,8 @@ __attribute__((always_inline)) INLINE static void rt_clean(
* 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, &grackle_rates);
+ local_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]);
diff --git a/src/rt/GEAR/rt_grackle_utils.h b/src/rt/GEAR/rt_grackle_utils.h
index 815e36ca8777bacd54d83c9b71721d3bba08633f..599221737365832fc1f8df665fe67757006a93e3 100644
--- a/src/rt/GEAR/rt_grackle_utils.h
+++ b/src/rt/GEAR/rt_grackle_utils.h
@@ -19,6 +19,9 @@
#ifndef SWIFT_RT_GRACKLE_UTILS_H
#define SWIFT_RT_GRACKLE_UTILS_H
+/* skip deprecation warnings. I cleaned old API calls. */
+#define OMIT_LEGACY_INTERNAL_GRACKLE_FUNC
+
/* need hydro gamma */
#include "hydro.h"
@@ -45,6 +48,7 @@
**/
__attribute__((always_inline)) INLINE static void rt_init_grackle(
code_units *grackle_units, chemistry_data *grackle_chemistry_data,
+ chemistry_data_storage *grackle_chemistry_rates,
float hydrogen_mass_fraction, const int grackle_verb,
const int case_B_recombination, const struct unit_system *us) {
@@ -66,23 +70,24 @@ __attribute__((always_inline)) INLINE static void rt_init_grackle(
/* Chemistry Parameters */
/* -------------------- */
- /* More details on https://grackle.readthedocs.io/en/latest/Parameters.html */
+ /* More details on
+ * https://grackle.readthedocs.io/en/grackle-3.2.0/Integration.html#chemistry-data
+ */
- /* TODO: cleanup later with all other grackle stuff */
- /* rtp->grackle_chemistry_data = _set_default_chemistry_parameters(); */
- if (set_default_chemistry_parameters(grackle_chemistry_data) == 0) {
+ if (local_initialize_chemistry_parameters(grackle_chemistry_data) == 0) {
error("Error in set_default_chemistry_parameters.");
}
+
/* chemistry on */
grackle_chemistry_data->use_grackle = 1;
/* cooling on */
/* NOTE: without cooling on, it also won't heat... */
grackle_chemistry_data->with_radiative_cooling = 1;
- /* 6 species atomic H and He-> */
+ /* 6 species atomic H and He */
grackle_chemistry_data->primordial_chemistry = 1;
- /* dust processes */
+ /* No dust processes */
grackle_chemistry_data->dust_chemistry = 0;
- /* H2 formation on dust */
+ /* No H2 formation on dust */
grackle_chemistry_data->h2_on_dust = 0;
/* metal cooling (uses Cloudy) off (for now) */
grackle_chemistry_data->metal_cooling = 0;
@@ -90,27 +95,21 @@ __attribute__((always_inline)) INLINE static void rt_init_grackle(
grackle_chemistry_data->cmb_temperature_floor = 1;
/* UV background off */
grackle_chemistry_data->UVbackground = 0;
- /* data file - currently not used-> */
+ /* data file - currently not used */
grackle_chemistry_data->grackle_data_file = "";
/* adiabatic index */
grackle_chemistry_data->Gamma = hydro_gamma;
/* we'll provide grackle with ionization and heating rates from RT */
grackle_chemistry_data->use_radiative_transfer = 1;
+
/* fraction by mass of Hydrogen in the metal-free portion of the gas */
grackle_chemistry_data->HydrogenFractionByMass = hydrogen_mass_fraction;
/* Use case B recombination? (On-the-spot approximation) */
grackle_chemistry_data->CaseBRecombination = case_B_recombination;
- /* TODO: cleanup later with all other grackle stuff */
- /* Initialize the chemistry_data_storage object to be
- * able to use local functions */
- /* chemistry_data_storage* chem_data_storage = */
- /* malloc(sizeof(chemistry_data_storage)); */
- /* rtp->grackle_chemistry_rates = chem_data_storage; */
- /* if (!_initialize_chemistry_data(&rtp->grackle_chemistry_data, */
- /* rtp->grackle_chemistry_rates, */
- /* &rtp->grackle_units)) { */
- if (initialize_chemistry_data(grackle_units) == 0) {
+ if (local_initialize_chemistry_data(grackle_chemistry_data,
+ grackle_chemistry_rates,
+ grackle_units) == 0) {
error("Error in initialize_chemistry_data");
}
}
@@ -177,6 +176,8 @@ rt_get_grackle_particle_fields(grackle_field_data *grackle_fields,
grackle_fields->HDI_density = NULL;
/* for metal_cooling = 1 */
grackle_fields->metal_density = NULL;
+ /* for use_dust_density_field = 1 */
+ grackle_fields->dust_density = NULL;
/* volumetric heating rate (provide in units [erg s^-1 cm^-3]) */
grackle_fields->volumetric_heating_rate = NULL;
@@ -196,6 +197,10 @@ rt_get_grackle_particle_fields(grackle_field_data *grackle_fields,
* (provide in units [erg s^-1 cm^-3] / nHI in cgs) */
grackle_fields->RT_heating_rate = malloc(FIELD_SIZE * sizeof(gr_float));
+ grackle_fields->H2_self_shielding_length = NULL;
+ grackle_fields->H2_custom_shielding_factor = NULL;
+ grackle_fields->isrf_habing = NULL;
+
for (int i = 0; i < FIELD_SIZE; i++) {
grackle_fields->density[i] = density;
@@ -249,24 +254,37 @@ __attribute__((always_inline)) INLINE static void rt_clean_grackle_fields(
free(grackle_fields->grid_start);
free(grackle_fields->grid_end);
+ /* initial quantities */
free(grackle_fields->density);
free(grackle_fields->internal_energy);
/* free(grackle_fields->x_velocity); */
/* free(grackle_fields->y_velocity); */
/* free(grackle_fields->z_velocity); */
+
+ /* for primordial_chemistry >= 1 */
free(grackle_fields->HI_density);
free(grackle_fields->HII_density);
free(grackle_fields->HeI_density);
free(grackle_fields->HeII_density);
free(grackle_fields->HeIII_density);
free(grackle_fields->e_density);
+
+ /* for primordial_chemistry >= 2 */
/* free(grackle_fields->HM_density); */
/* free(grackle_fields->H2I_density); */
/* free(grackle_fields->H2II_density); */
+
+ /* for primordial_chemistry >= 3 */
/* free(grackle_fields->DI_density); */
/* free(grackle_fields->DII_density); */
/* free(grackle_fields->HDI_density); */
+
+ /* for metal_cooling = 1 */
/* free(grackle_fields->metal_density); */
+
+ /* for use_dust_density_field = 1 */
+ /* free(grackle_fields->dust_density); */
+
/* free(grackle_fields->volumetric_heating_rate); */
/* free(grackle_fields->specific_heating_rate); */
@@ -275,6 +293,10 @@ __attribute__((always_inline)) INLINE static void rt_clean_grackle_fields(
free(grackle_fields->RT_HeII_ionization_rate);
free(grackle_fields->RT_H2_dissociation_rate);
free(grackle_fields->RT_heating_rate);
+
+ /* free(grackle_fields->H2_self_shielding_length); */
+ /* free(grackle_fields->H2_custom_shielding_factor); */
+ /* free(grackle_fields->isrf_habing); */
}
/**
diff --git a/src/rt/GEAR/rt_properties.h b/src/rt/GEAR/rt_properties.h
index a9226cc7568a0458d8cbe8998c9c2a216c333d13..49b0c9c60d50067f7927faa22bcc6df5f549db2d 100644
--- a/src/rt/GEAR/rt_properties.h
+++ b/src/rt/GEAR/rt_properties.h
@@ -73,7 +73,7 @@ struct rt_props {
float mass_fraction_HeI_init;
float mass_fraction_HeII_init;
float mass_fraction_HeIII_init;
- float number_density_electrons_init; /* todo: do we need this? */
+ /* float number_density_electrons_init; [> todo: do we need this? <] */
/* Hydrogen and Helium mass fractions of the non-metal portion of the gas */
float hydrogen_mass_fraction;
@@ -114,17 +114,16 @@ struct rt_props {
/*! grackle chemistry data */
chemistry_data grackle_chemistry_data;
- /* use case B recombination? */
+ /*! grackle chemistry data storage
+ * (needed for local function calls) */
+ chemistry_data_storage grackle_chemistry_rates;
+
+ /*! use case B recombination? */
int case_B_recombination;
- /* make grackle talkative? */
+ /*! make grackle talkative? */
int grackle_verbose;
- /* TODO: cleanup later with all other grackle stuff */
- /*! grackle chemistry data storage
- * (needed for local function calls) */
- /* chemistry_data_storage* grackle_chemistry_rates; */
-
#ifdef SWIFT_RT_DEBUG_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
@@ -452,8 +451,8 @@ __attribute__((always_inline)) INLINE static void rt_props_init(
rtp->case_B_recombination = parser_get_opt_param_int(
params, "GEARRT:case_B_recombination", /*default=*/1);
rt_init_grackle(&rtp->grackle_units, &rtp->grackle_chemistry_data,
- rtp->hydrogen_mass_fraction, rtp->grackle_verbose,
- rtp->case_B_recombination, us);
+ &rtp->grackle_chemistry_rates, rtp->hydrogen_mass_fraction,
+ rtp->grackle_verbose, rtp->case_B_recombination, us);
/* Pre-compute interaction rates/cross sections */
/* -------------------------------------------- */
@@ -501,6 +500,7 @@ __attribute__((always_inline)) INLINE static void rt_struct_restore(
"RT properties struct");
/* Set up stuff that needs array allocation */
rt_init_grackle(&props->grackle_units, &props->grackle_chemistry_data,
+ &props->grackle_chemistry_rates,
props->hydrogen_mass_fraction, props->grackle_verbose,
props->case_B_recombination, us);
diff --git a/src/rt/GEAR/rt_thermochemistry.h b/src/rt/GEAR/rt_thermochemistry.h
index 917016a48adca7a29e254549daa03697721b6e4e..273a240e3d7ed0ec89f6e6928224ff5d25451409 100644
--- a/src/rt/GEAR/rt_thermochemistry.h
+++ b/src/rt/GEAR/rt_thermochemistry.h
@@ -65,8 +65,12 @@ __attribute__((always_inline)) INLINE static void rt_tchem_first_init_part(
p->rt_data.tchem.mass_fraction_HeIII = rt_props->mass_fraction_HeIII_init;
}
+ /* pretend you have nonzero density so the check doesn't reset the mass
+ * fractions */
+ p->rho = 1.f;
/* Check that we didn't do something stupid */
rt_check_unphysical_mass_fractions(p);
+ p->rho = 0.f;
/* Check that the Hydrogen and Helium mass fractions correspond to those
* provided by the user in the parameter file. This mass fraction is also
@@ -151,9 +155,9 @@ INLINE static void rt_do_thermochemistry(
/* Note: `grackle_rates` is a global variable defined by grackle itself.
* Using a manually allocd and initialized variable here fails with MPI
* for some reason. */
- if (local_solve_chemistry(&rt_props->grackle_chemistry_data, &grackle_rates,
- &rt_props->grackle_units, &particle_grackle_data,
- dt) == 0)
+ if (local_solve_chemistry(
+ &rt_props->grackle_chemistry_data, &rt_props->grackle_chemistry_rates,
+ &rt_props->grackle_units, &particle_grackle_data, dt) == 0)
error("Error in solve_chemistry.");
/* copy updated grackle data to particle */
@@ -282,9 +286,9 @@ __attribute__((always_inline)) INLINE static float rt_tchem_get_tchem_time(
* Using a manually allocd and initialized variable here fails with MPI
* for some reason. */
gr_float tchem_time;
- if (local_calculate_cooling_time(&rt_props->grackle_chemistry_data,
- &grackle_rates, &rt_props->grackle_units,
- &particle_grackle_data, &tchem_time) == 0)
+ if (local_calculate_cooling_time(
+ &rt_props->grackle_chemistry_data, &rt_props->grackle_chemistry_rates,
+ &rt_props->grackle_units, &particle_grackle_data, &tchem_time) == 0)
error("Error in calculate_cooling_time.");
/* Clean up after yourself. */
diff --git a/src/rt/GEAR/rt_unphysical.h b/src/rt/GEAR/rt_unphysical.h
index 83c31063d06bf95c0b009e53078b088db4b5f23d..39bb0086d4b77697bb5dd7959e4c7d9f2652231e 100644
--- a/src/rt/GEAR/rt_unphysical.h
+++ b/src/rt/GEAR/rt_unphysical.h
@@ -183,13 +183,14 @@ rt_check_unphysical_mass_fractions(struct part* restrict p) {
* particles, while they themselves remain inactive. The density of such
* inactive particles however remains zero until the particle is active
* again. See issue #833. */
- if (p->conserved.mass * p->rho <= 0.f) {
+
+ if (p->conserved.mass <= 0.f || p->rho <= 0.f) {
/* Deal with unphysical situations and vacuum. */
- p->rt_data.tchem.mass_fraction_HI = 0.f;
- p->rt_data.tchem.mass_fraction_HII = 0.f;
- p->rt_data.tchem.mass_fraction_HeI = 0.f;
- p->rt_data.tchem.mass_fraction_HeII = 0.f;
- p->rt_data.tchem.mass_fraction_HeIII = 0.f;
+ p->rt_data.tchem.mass_fraction_HI = RT_GEAR_TINY_MASS_FRACTION;
+ p->rt_data.tchem.mass_fraction_HII = RT_GEAR_TINY_MASS_FRACTION;
+ p->rt_data.tchem.mass_fraction_HeI = RT_GEAR_TINY_MASS_FRACTION;
+ p->rt_data.tchem.mass_fraction_HeII = RT_GEAR_TINY_MASS_FRACTION;
+ p->rt_data.tchem.mass_fraction_HeIII = RT_GEAR_TINY_MASS_FRACTION;
return;
}
diff --git a/src/runner_ghost.c b/src/runner_ghost.c
index 48c8055b8cfb8287d649b0036da4af67c01eba85..091589f52919b7a370ceaae5809fb16ef47c066e 100644
--- a/src/runner_ghost.c
+++ b/src/runner_ghost.c
@@ -327,7 +327,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
if ((h_new == left[i] && h_old == right[i]) ||
(h_old == left[i] && h_new == right[i])) {
- /* Bissect the remaining interval */
+ /* Bisect the remaining interval */
sp->h = pow_inv_dimension(
0.5f * (pow_dimension(left[i]) + pow_dimension(right[i])));
@@ -748,7 +748,7 @@ void runner_do_black_holes_density_ghost(struct runner *r, struct cell *c,
if ((h_new == left[i] && h_old == right[i]) ||
(h_old == left[i] && h_new == right[i])) {
- /* Bissect the remaining interval */
+ /* Bisect the remaining interval */
bp->h = pow_inv_dimension(
0.5f * (pow_dimension(left[i]) + pow_dimension(right[i])));
@@ -1357,7 +1357,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
if ((h_new == left[i] && h_old == right[i]) ||
(h_old == left[i] && h_new == right[i])) {
- /* Bissect the remaining interval */
+ /* Bisect the remaining interval */
p->h = pow_inv_dimension(
0.5f * (pow_dimension(left[i]) + pow_dimension(right[i])));
diff --git a/src/runner_time_integration.c b/src/runner_time_integration.c
index 80ba9e2912339d367127ed450cd2f32442fc18e1..152fc8ae342e6b2bdf27b486443e0eb0f07cf152 100644
--- a/src/runner_time_integration.c
+++ b/src/runner_time_integration.c
@@ -1752,7 +1752,12 @@ void runner_do_collect_rt_times(struct runner *r, struct cell *c,
* rt_advanced_cell_time should be called exactly once before
* collect times. Except on the first subcycle, because the
* collect_rt_times task shouldn't be called in the main steps.
- * In that case, it should be exactly 2. */
+ * In that case, it should be exactly 2.
+ * This is only valid if the cell has been active this step.
+ * Otherwise, rt_advance_cell_time will not have run, yet the
+ * rt_collect_cell_times call may still be executed if the top
+ * level cell is above the super level cell. */
+ if (!cell_is_rt_active(c, e)) return;
if (e->ti_current_subcycle - c->rt.ti_rt_end_min == e->ti_current) {
/* This is the first subcycle */
if (c->rt.advanced_time != 2)
diff --git a/src/statistics.c b/src/statistics.c
index 87212bda610a9e3cec1b12b6e058054271fdd951..18bf27bb43e1a4b668a215a1a99c2af10bfcd0f9 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -183,9 +183,8 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
/* Collect metal mass */
stats.gas_Z_mass += chemistry_get_total_metal_mass_for_stats(p);
-#if defined(CHEMISTRY_EAGLE) || defined(CHEMISTRY_PS2020)
-#if defined(COOLING_EAGLE) || defined(COOLING_PS2020) || \
- defined(COOLING_CHIMES) || defined(COOLING_CHIMES_HYBRID)
+#if defined(CHEMISTRY_EAGLE)
+#if defined(COOLING_EAGLE) || defined(COOLING_PS2020)
const struct unit_system *us = e->internal_units;
const struct hydro_props *hydro_props = e->hydro_properties;
diff --git a/swift.c b/swift.c
index 4201c7117a23ece1d85649844486f97c9021d872..f877b0b16fac977610443aacbbe29126fb5140c4 100644
--- a/swift.c
+++ b/swift.c
@@ -644,6 +644,9 @@ int main(int argc, char *argv[]) {
error("Running with radiative transfer but compiled without it!");
}
#else
+ if (!with_rt) {
+ error("Running without radiative transfer but compiled with it!");
+ }
if (with_rt && !with_hydro) {
error(
"Error: Cannot use radiative transfer without gas, --hydro must be "
diff --git a/theory/RadiativeTransfer/GEARRT/1-meshless.tex b/theory/RadiativeTransfer/GEARRT/1-meshless.tex
new file mode 100644
index 0000000000000000000000000000000000000000..f1e3b95d2f1baef6ac0bef6980038c58cca40cee
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/1-meshless.tex
@@ -0,0 +1,115 @@
+%===================================================================================
+\section{Finite Volume Particle Method: Essentials} \label{chap:meshless}
+%===================================================================================
+
+The Finite Volume Particle methods (FVPM) are a way to solve hyperbolic systems of conservation
+laws of the form
+
+\begin{equation}
+ \DELDT{\mathcal{U}} + \nabla \mathcal{F}(\mathcal{U}) = 0 \label{eq:conservation-law}
+\end{equation}
+
+for some state vector $\mathcal{U}$ and flux tensor $\mathcal{F}(\mathcal{U})$
+and are used to solve the hydrodynamics equations. Since the M1 RT method can be written in the same
+form (see Section~\ref{chap:rt-equations}), the method will also be used to solve these equations.
+For completeness and quick reference, the most important formulae are listed below. For the actual
+implementation details, please refer to the `Gizmo' theory section in the SWIFT repository, and
+\citet{hopkinsGIZMONewClass2015}. For a step-by-step derivation and introduction, see
+\citet{ivkovicThesis}.
+
+
+
+The FVPM schemes rely on the ``partitions of unity''. The partition of unity assigned to a particle
+$i$ at some point $\x$ is defined as:
+
+\begin{align}
+ \psi_i(\x) &= \frac{1}{\omega(\x)} W(\x - \x_i, h(\x)) \label{eq:psi} \\
+ \omega(\x) &= \sum_j W(\x - \x_j, h(\x)) \label{eq:omega}
+\end{align}
+
+where $h(\x)$ is the smoothing length at $\x$ and $\omega(\x)$ is used to normalise the volume
+partition at any point $\x$.
+
+The associated particle volumes are given by
+
+\begin{align}
+ V_i &= \int_V \psi_i(\x) \de V \approx \frac{1}{\omega(\x_i)} \label{eq:volume}
+\end{align}
+
+
+The actual state updates are given by
+
+\begin{equation}
+ \frac{\D}{\D t} (V_i \mathcal{U}_{k,i}) + \sum_j \mathcal{F}_{k,ij} \cdot \Aijm = 0
+\end{equation}
+
+
+with
+
+\begin{equation}
+ \Aijm^\alpha = V_i \psitilde_j^\alpha (\x_i) - V_j\psitilde_i^\alpha (\x_j) \label{eq:HopkinsAij}
+\end{equation}
+
+for every component $k$ of the state vector $\mathcal{U}$ and for every coordinate dimension
+$\alpha$.
+
+The $\psitilde(\x)$ come from the $\mathcal{O}(h^2)$ accurate discrete gradient expression from
+\cite{lansonRenormalizedMeshfreeSchemes2008}:
+
+\begin{align}
+ \frac{\del}{\del x_{\alpha}} f(\x) \big{|}_{\x_i} &=
+ \sum_j \left( f(\x_j) - f(\x_i) \right) \psitilde_j^\alpha (\x_i) \label{eq:gradient} \\
+ \psitilde_j^\alpha (\x_i) &= \sum_{\beta = 1}^{\beta=\nu} \mathcal{B}_i^{\alpha \beta}
+ (\x_j - \x_i)^\beta \psi_j(\x_i) \\
+ \mathcal{B}_i &= \mathcal{E_i} ^ {-1} \label{eq:matrix-B}\\
+ \mathcal{E}_i^{\alpha \beta} &= \sum_j (\x_j - \x_i)^\alpha (\x_j - \x_i)^\beta \psi_j(\x_i)
+\label{eq:matrix-E}
+\end{align}
+
+
+where $\alpha$ and $\beta$ again represent the coordinate components for $\nu$ dimensions.
+
+
+The final explicit update formula is given by
+
+\begin{align}
+ \U_i^{n+1} = \U_i^{n} + \frac{\Delta t}{V_i} \sum_j \F_{\alpha, ij} \Aijm^\alpha
+\label{eq:meshless-explicit}
+\end{align}
+
+
+
+
+
+
+
+
+%===================================================================================
+\subsection{Individual Timestepping}
+%===================================================================================
+
+With particles being allowed to have individual time steps, the update formula for the finite volume
+particle methods (eq.~\ref{eq:meshless-explicit}) needs to be slightly updated in order to
+maintain its conservative properties. In particular, rather than exchanging fluxes, particles now
+need to exchange time integrated fluxes, which are integrated using a trivial explicit Euler scheme:
+
+\begin{align}
+\U_i^{n+1} =
+ \U_i^n + \frac{1}{V_i} \sum_j \left[
+ \min\{ \Delta t_i, \Delta t_j \} \F_{\alpha, ij} \Aijm^\alpha
+ \right]
+\label{eq:meshless-explicit-individual-timesteps}
+\end{align}
+
+The time integration needs to be the smaller of the two time step sizes $\{\Delta t_i,\ \Delta
+t_j\}$ since if one particle has a smaller time step size than the other, it also means that there
+will be several interactions between that particle pair, each with the time step size of the smaller
+of the two. Suppose for example $\Delta t_i = 4 \Delta t_j$. Then for a single update of particle
+$i$, there will be in total 4 updates of particle $j$. Each update of particle $j$ will entail a
+symmetric (time integrated) flux exchange with particle $i$. Therefore the time integration of the
+flux to be exchanged must be applied over the smaller time interval $\Delta t_j$ for each of
+the four (time integrated) flux exchanges.
+
+
+
+
diff --git a/theory/RadiativeTransfer/GEARRT/2-rt-equations.tex b/theory/RadiativeTransfer/GEARRT/2-rt-equations.tex
new file mode 100644
index 0000000000000000000000000000000000000000..ff8e1335ad4083c777185d4d853fedf733fac0b2
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/2-rt-equations.tex
@@ -0,0 +1,224 @@
+%=============================================================================
+\section{The Equations of Moment-Based Radiative Transfer}\label{chap:rt-equations}
+%=============================================================================
+
+
+%-----------------------------------------------------------------------
+\subsection{The Equations of Radiative Transfer and the M1 Closure}
+%-----------------------------------------------------------------------
+
+
+%-----------------------------------------------------------------------
+\subsubsection{The Equations of Radiative Transfer}
+%-----------------------------------------------------------------------
+
+
+Radiative transfer (RT) and radiation hydrodynamics (RHD) contain a plethora of variables and
+coefficients. For clarity, an overview of the relevant quantities and coefficients is given
+in Appendix~\ref{app:variables} along with their respective units.
+
+
+The equation of radiative transfer is given by:
+
+\begin{align}
+ \frac{1}{c} \DELDT{I_\nu} + \mathbf{n} \cdot \nabla I_\nu
+ &= \eta_\nu - \alpha_\nu I_\nu \nonumber \\
+ &= \eta_\nu - \sum_j^{\text{photo-absorbing\ species}} \sigma_{j,\nu} n_j I_\nu \ .
+\label{eq:RT-sigma}
+\end{align}
+
+\begin{itemize}
+\item $I_\nu$ is the specific intensity and has units of erg cm$^{-3}$ rad$^{-1}$ Hz$^{-1}$
+s$^{-1}$.
+\item $\eta_\nu$ is a source function of radiation, i.e. the term describing radiation being added
+along the (dimensionless) direction $\mathbf{n}$ due to yet unspecified processes, and has units of
+erg cm$^{-3}$ rad$^{-1}$ Hz$^{-1}$ s$^{-1}$, which is the same as the units of the specific
+intensity $I_\nu$ per cm.
+\item $\alpha_\nu$ is an absorption coefficient, describes how much radiation is being removed, and
+has units of cm$^{-1}$. Naturally only as much radiation as is currently present can be removed, and
+so the sink term must be proportional to the local specific intensity $I_\nu$.
+\end{itemize}
+
+The equation holds for any photon frequency $\nu$ individually. In eq.~\ref{eq:RT-sigma} we split
+the absorption coefficient $\alpha_\nu$ into the sum over the photo-absorbing species $j$, which
+for GEAR-RT will only be the main constituents of primordial gas, namely hydrogen, helium,
+and singly ionized helium. The photo-absorption process is expressed via interaction cross sections
+$\sigma_{j,\nu}$, which has units of cm$^2$, while $n_j$ represents the number density of
+photo-absorbing species $j$ in cm$^{-3}$.
+
+
+
+%-----------------------------------------------------------------------
+\subsubsection{Moments of The Equations of Radiative Transfer}
+%-----------------------------------------------------------------------
+
+
+GEAR-RT solves for the (angular) moments of the equation of radiative transfer, which are obtained
+through integrating eq.~\ref{eq:RT-sigma} over the entire solid angle to obtain the zeroth moment
+equation, and over the entire solid angle multiplied by the direction unit vector $\mathbf{n}$ for
+the first moment equation.
+Additionally, we make use of the following quantities:
+
+\begin{align*}
+ E_\nu (\x, t) &= \int_{4 \pi} \frac{I_\nu}{c} \de \Omega
+ && \text{total energy density }
+ && [E_\nu] = \frac{\text{erg}}{\text{cm}^3 \text{ Hz}}\\
+ \Fbf_\nu(\x, t) &= \int_{4 \pi} I_\nu \mathbf{n} \de \Omega
+ && \text{radiation flux }
+ && [\Fbf_\nu] = \frac{\text{erg}}{\text{cm}^2 \text{ s Hz}}\\
+ \mathds{P}_\nu (\x, t) &= \int_{4 \pi} \frac{I_\nu}{c} \mathbf{n} \otimes \mathbf{n} \de \Omega
+ && \text{radiation pressure tensor }
+ && [\mathds{P}_\nu ] = \frac{\text{erg}}{\text{cm}^3 \text{ Hz}}
+\end{align*}
+
+where $\mathbf{n} \otimes \mathbf{n}$ denotes the outer product, which in components $k$, $l$ gives
+%
+\begin{align*}
+ (\mathbf{n} \otimes \mathbf{n})_{kl} = \mathbf{n}_k \mathbf{n}_l \ .
+\end{align*}
+
+
+
+This gives us the following equations:
+
+
+\begin{align}
+ \DELDT{E_\nu} + \nabla \cdot \Fbf_\nu &=
+ - \sum\limits_{j}^{\absorbers} n_j \sigma_{\nu j} c E_\nu + \dot{E}_\nu
+ \label{eq:dEdt-freq} \\
+ \DELDT{\Fbf_\nu} + c^2 \ \nabla \cdot \mathds{P}_\nu &=
+ - \sum\limits_{j}^{\absorbers} n_j \sigma_{\nu j} c \Fbf_\nu
+ \label{eq:dFdt-freq}
+\end{align}
+
+
+
+Note that $E_\nu$ (and $\dot{E}_\nu$) is the radiation energy \emph{density} (and \emph{density}
+injection rate) in the frequency interval between frequency $\nu$ and $\nu + \de \nu$ and has
+units of $\text{erg / cm}^3 \text{ / Hz}$ (and $\text{erg / cm}^3 \text{ / Hz / s}$).
+$\Fbf$ is the radiation flux, and has units of $\text{erg / cm}^2 \text{ / Hz / s}$, i.e.
+dimensions of energy per area per frequency per time.
+
+Furthermore, it is assumed that the source term $\dot{E}_\nu$ stems from point sources which radiate
+isotropically. This assumption has the consequence that the vector net flux $\Fbf_\nu$ must sum up
+to zero, and hence the corresponding source terms in eq.~\ref{eq:dFdt-freq} are zero.
+
+
+
+%-----------------------------------------------------------------------
+\subsubsection{The M1 Closure}
+%-----------------------------------------------------------------------
+
+
+
+To close this set of equations, a model for the pressure tensor $\mathds{P}_\nu$ is necessary. We
+use the so-called ``M1 closure'' \citep{levermoreRelatingEddingtonFactors1984a} where we describe
+the pressure tensor via the Eddington tensor $\mathds{D}_\nu$:
+
+
+\begin{equation*}
+ \mathds{P}_\nu = \mathds{D}_\nu E_\nu
+\end{equation*}
+
+The Eddington tensor is a dimensionless quantity that encapsulates the local radiation field
+geometry and its effect in the radiation flux conservation equation. The M1 closure sets the
+Eddington tensor to have the form:
+
+\begin{align*}
+\mathds{D}_\nu &=
+ \frac{1- \chi_\nu}{2} \mathds{I} + \frac{3 \chi_\nu - 1}{2} \mathbf{n}_\nu \otimes
+ \mathbf{n}_\nu
+ \label{eq:eddington-freq} \\
+\mathbf{n}_\nu &=
+ \frac{\Fbf_\nu}{|\Fbf_\nu|} \\
+\chi_\nu &=
+ \frac{3 + 4 f_\nu ^2}{5 + 2 \sqrt{4 - 3 f_\nu^2}} \\
+f_\nu &=
+ \frac{|\Fbf_\nu|}{c E_\nu}
+\end{align*}
+
+
+
+
+
+
+
+
+
+
+
+
+%---------------------------------------------------------------
+\subsection{Photo-ionization and Photo-heating Rates}
+%---------------------------------------------------------------
+
+
+In the context of radiative transfer and photo-ionization, the photo-ionization rate $\Gamma_{\nu,
+j}$ in units of s$^{-1}$ for photons with frequency $\nu$ and a photo-ionizing particle species $j$
+is then given by
+
+\begin{align}
+ \DELDT{n_j} = -\Gamma_{\nu, j} \ n_j = - c \ \sigma_{\nu j} \ N_\nu \ n_j
+\end{align}
+
+where $N_\nu = E_\nu / (h \nu)$ is the photon number density, and $n_j$ is the number density of
+the photo-ionizing particle species. Note that both the interaction cross sections and the
+photo-ionizing species $j$ are specific to a frequency $\nu$.
+
+For the cross sections, we use the analytic fits for the photo-ionization cross sections from
+\cite{vernerAtomicDataAstrophysics1996} (via \cite{ramses-rt13}), which are given by
+
+\begin{align}
+\sigma(E) &= \sigma_0 F(y) \times 10^{-18} \text{ cm}^2 \label{eq:sigma-parametrizaiton}
+\\
+F(y) &= \left[(x - 1)^2 + y_w^2 \right] y ^{0.5 P - 5.5} \left( 1 + \sqrt{y / y_a} \right)^{-P}
+\\
+x &= \frac{E}{E_0} - y_0 \\
+y &= \sqrt{x^2 + y_1^2}
+\end{align}
+
+where $E$ is the photon energy $E = h \nu$ in eV, and $\sigma_0$, $E_0$, $y_w$, $y_a$, $P$, $y_0$,
+and $y_1$ are fitting parameters. The fitting parameter values for hydrogen, helium, and singly
+ionized helium are given in Table~\ref{tab:cross-sections}. The thresholds are given as frequencies
+in eqs.~\ref{eq:nuIonHI}-\ref{eq:nuIonHeII}. Below this threshold, no ionization can take place, and
+hence the cross sections are zero.
+
+
+\input{tables/fitting_parameters}
+
+
+
+
+
+
+Conversely, the rate at which photons are absorbed, i.e. ``destroyed'', must be equal to the
+photo-ionization rate, which means
+
+\begin{align*}
+\DELDT{E_\nu} &= h \ \nu \DELDT{N_\nu} = -h\ \nu \ c \ \sigma_{\nu j} n_j N_\nu
+\end{align*}
+
+
+
+Finally, the photo-heating rate is modeled as the rate of excess energy absorbed by the gas during
+photo-ionizing collisions. To ionize an atom, the photons must carry a minimal energy corresponding
+to the ionizing frequency $\nu_{ion,j}$ for a photo-ionizing species $j$. In the case of hydrogen
+and helium, their values are
+
+\begin{align}
+ \nu_{\text{ion,HI}} &= 2.179 \times 10^{-11} \text{ erg} = 13.60 \text{ eV} \label{eq:nuIonHI}\\
+ \nu_{\text{ion,HeI}} &= 3.940 \times 10^{-11} \text{ erg} = 24.59 \text{ eV}
+\label{eq:nuIonHeI}\\
+ \nu_{\text{ion,HeII}} &= 8.719 \times 10^{-11} \text{ erg} = 54.42 \text{ eV}
+\label{eq:nuIonHeII}
+\end{align}
+
+All excess energy is added to the gas in the form of internal energy, and the heating rate
+$\mathcal{H}$ (in units of erg cm$^{-3}$ s$^{-1}$ ) for photons of some frequency $\nu$ and some
+photo-ionizing species $j$ is described by
+
+\begin{align*}
+\mathcal{H}_{\nu, j} = (h \nu - h \nu_{ion,j}) \ c \ \sigma_{\nu j} \ n_j \ N_\nu
+\end{align*}
+
+
diff --git a/theory/RadiativeTransfer/GEARRT/3-rt-numerical-method.tex b/theory/RadiativeTransfer/GEARRT/3-rt-numerical-method.tex
new file mode 100644
index 0000000000000000000000000000000000000000..a3350bc7a62bb00ada5d9c0e7922304bfc094274
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/3-rt-numerical-method.tex
@@ -0,0 +1,857 @@
+%======================================================
+\section{Numerical Solution Strategy}\label{chap:rt-numerical-strategy}
+%======================================================
+
+
+%---------------------------------------------
+\subsection{Discretization of Frequencies}
+%---------------------------------------------
+
+
+The moments of the equation of radiative transfer, given in
+eqns.~\ref{eq:dEdt-freq}~-~\ref{eq:dFdt-freq}, are still frequency specific, and would need to be
+solved for each frequency between 0 Hz and infinity individually. This is obviously not a feasible
+approach. Instead, these equations need to be discretized in frequency first. I follow the approach
+of \cite{ramses-rt13} and for a rough approximation of multi-frequency, a relevant frequency range
+is split into a number $M$ of mutually exclusive groups of frequency ranges:
+
+\begin{align*}
+ & [\nu_{00}, \nu_{01}\ : \ \nu_{10}, \nu_{11}\ : ... \ : \ \nu_{M0}, \nu_{M1}\ ] =
+ [\nu_{0}, \infty [
+\end{align*}
+
+The frequency ranges are chosen to be convenient for us. Specifically, since the interaction cross
+sections of ionizing species are zero below the ionizing frequency of their corresponding species,
+it makes sense to use the various ionizing frequencies given in
+eqs.~\ref{eq:nuIonHI}-\ref{eq:nuIonHeII} as the boundaries for the frequency groups.
+
+
+
+
+Rather than treating the photon energy densities and fluxes for individual frequencies, we now use
+their integrated averages over a frequency group. For any frequency interval, or group, $i$, we have
+
+\begin{align}
+ E_i &=
+ \int\limits_{\nu_{i0}}^{\nu_{i1}} E_\nu d\nu \\
+ \Fbf_i &=
+ \int\limits_{\nu_{i0}}^{\nu_{i1}} \Fbf_\nu d\nu
+\end{align}
+
+
+giving us the following equations (discretized in frequency) to solve:
+
+
+\begin{align}
+ \DELDT{E_i} + \nabla \cdot \Fbf_i &=
+ - \sum\limits_{j}^{\absorbers} n_j \sigma_{i j}^N c E_i + \dot{E}_i^* + \dot{E}_i ^{rec}
+ \label{eq:dEdt-group} \\
+ \DELDT{\Fbf_i} + c^2 \ \nabla \cdot \mathds{P}_i &=
+ - \sum\limits_{j}^{\absorbers} n_j \sigma_{i j}^N c \Fbf_i
+ \label{eq:dFdt-group}
+\end{align}
+
+
+The expression for the number weighted average cross section $\sigma_{ij}^N$ is given in eq.
+\ref{eq:sigma_N}. The radiation pressure tensor is discretized in the same manner:
+
+\begin{align}
+ \mathds{P}_i &=
+ \mathds{D}_i E_i \label{eq:pressure-tensor-group-start}\\
+ \mathds{D}_i &=
+ \frac{1- \chi_i}{2} \mathds{I} + \frac{3 \chi_i - 1}{2} \mathbf{n}_i \otimes \mathbf{n}_i
+ \label{eq:eddington-group} \\
+ \mathbf{n}_i &=
+ \frac{\Fbf_i}{|\Fbf_i|} \\
+ \chi_i &=
+ \frac{3 + 4 f_i ^2}{5 + 2 \sqrt{4 - 3 f_i^2}} \\
+ f_i &=
+ \frac{|\Fbf_i|}{c E_i} \label{eq:pressure-tensor-group-end}
+\end{align}
+
+
+
+For the computation of the photo-heating and photo-ionization rates, which will be discussed
+subsequently, we need to introduce the mean photon energy $\overline{\epsilon}_i$ of the frequency
+bin $i$. The exact value cannot be known, as it would require the knowledge of the energy density at
+each frequency, which is not available as we are using quantities averaged over frequency intervals.
+Instead, we guess a reasonable spectrum. \cite{ramses-rt13} recommend taking an average value among
+all stellar radiation sources. GEAR-RT currently only works with a blackbody spectrum, which for a
+characteristic temperature $T_{bb}$ is given by:
+
+\begin{align}
+ \mathcal{J}_\nu(T_{bb}) = \frac{2 \nu^2}{c^2} \frac{h \nu}{\exp\left(h\nu/k_B T_{bb}\right) -
+1} \ .
+ \label{eq:blackbody}
+\end{align}
+
+The temperature $T_{bb}$ in this case would be some effective temperature for stellar sources, not
+the local temperature of the gas. With an assumed spectral shape, the mean photon energy
+$\overline{\epsilon}_i$ in each frequency group $i$ can be estimated as
+
+
+\begin{align}
+\overline{\epsilon}_i \equiv
+ \frac{E_i}{N_i} =
+ \frac{
+ \int\limits_{\nu_{i0}}^{\nu_{i1}} E_\nu \ d\nu
+ }{
+ \int\limits_{\nu_{i0}}^{\nu_{i1}} N_\nu \ d\nu
+ }
+ \approx
+ \frac{
+ \int\limits_{\nu_{i0}}^{\nu_{i1}} \mathcal{J}_\nu \ d\nu
+ }{
+ \int\limits_{\nu_{i0}}^{\nu_{i1}} \mathcal{J}_\nu / h\nu \ d\nu
+ }
+\end{align}
+
+These integrals are performed using GSL integrator functions.
+
+
+Using the mean photon energy $\overline{\epsilon}_i$, the photo-heating rate $\mathcal{H}_{i,j}$
+of the gas for the photon group $i$ and ionizing species $j$ then becomes:
+
+\begin{align}
+\mathcal{H}_{i, j} &=
+\left[
+ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu h \nu N_\nu \sigma_{j\nu} -
+ h \nu_{ion,j}\
+ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu N_\nu \sigma_{j\nu} \
+\right] \ c \ n_j \\
+%
+&=
+\left[
+ \sigma_{ij}^E - h \nu_{ion,j}\ \sigma_{ij}^N /\ \overline{\epsilon}_i
+\right] E_i\ c \ n_j
+\label{eq:photoheating-group}
+\end{align}
+
+
+
+
+And the photo-ionization rate can be written as
+
+\begin{align}
+\Gamma_{i, j}
+ &=
+ c \ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ \sigma_{\nu j} N_\nu \\
+ &= c \ \sigma_{ij}^N N_i
+ = c \ \sigma_{ij}^N E_i / \overline{\epsilon}_i \label{eq:photoionization-group}
+\end{align}
+
+
+Finally, the photon absorption (destruction) rates for a frequency bin $i$ and ionizing species $j$
+are given by
+
+\begin{align}
+\DELDT{E_i} &=
+ \deldt{(N_i \ \overline{\epsilon}_i)} =
+ \overline{\epsilon}_i \DELDT{N_i} =
+ -\overline{\epsilon}_i c \ \sigma_{i j} ^ N n_j
+\end{align}
+
+
+Here we have introduced the number- and energy-weighted average cross sections:
+\begin{align}
+\sigma_{ij}^N &=
+ \frac{
+ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ N_\nu \ \sigma_{j\nu}
+ } {
+ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ N_\nu
+ }
+ \approx
+ \frac{
+ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ \mathcal{J}(\nu) / (h \ \nu) \sigma_{j\nu}
+ } {
+ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ \mathcal{J}(\nu) / (h \ \nu)
+ } \label{eq:sigma_N}\\
+\sigma_{ij}^E &=
+ \frac{
+ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ h \nu N_\nu \ \sigma_{j\nu}
+ } {
+ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ h \nu N_\nu
+ }
+ \approx
+ \frac{
+ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ \mathcal{J}(\nu) \ \sigma_{j\nu}
+ } {
+ \int\limits_{\nu_{i0}}^{\nu_{i1}}\de \nu \ \mathcal{J}(\nu)
+ } \label{eq:sigma_E}
+\end{align}
+
+
+
+
+
+
+
+
+%---------------------------------------------
+\subsection{One RT Time Step}
+%---------------------------------------------
+
+
+
+%---------------------------------------------
+\subsubsection{Outline}\label{chap:rt-numerics-outline}
+%---------------------------------------------
+
+
+For each photon frequency group, the equations \ref{eq:dEdt-group} and \ref{eq:dFdt-group} are
+solved with an operator-splitting strategy. Following the approach of \cite{ramses-rt13}, the
+equations are decomposed into three steps that are executed in sequence over the same time step
+$\Delta t$:
+
+\begin{enumerate}
+
+\item \textbf{Photon injection step}: the radiation from radiative sources is injected into the
+grid.
+
+\item \textbf{Photon transport step}: Photons are transported in space by solving the homogenized
+equations \ref{eq:dEdt-group} and \ref{eq:dFdt-group}, i.e. the right hand side of these equations
+is set to zero.
+
+\item \textbf{Thermochemistry step}: The rest of the source terms (the right hand side) of the
+equations \ref{eq:dEdt-group} and \ref{eq:dFdt-group} is solved.
+\end{enumerate}
+
+In what follows, the numerical and algorithmic aspects of these three steps are discussed
+individually. But first, let's list the fundamental approach upon which GEAR-RT is built:
+
+\begin{itemize}
+\item The underlying method to solve the hyperbolic conservation laws for hydrodynamics and for
+radiative transfer is a finite volume particle method.
+\item The particles are co-moving with the gas, and the motion of the particles is determined by
+the hydrodynamics.
+\item Radiation fields at particle positions are treated in a static frame of reference. (In
+contrast to the hydrodynamics, which are treated in a Lagrangian manner.)
+\item Particles are given individual time step sizes, for both the hydrodynamics and the radiative
+transfer individually.
+\item In a simulation step, the hydrodynamics are solved before the radiative transfer. Since RT is
+solved in a static frame of reference, this means that we can re-use all neighbour dependent data,
+in particular the smoothing lengths, partitions of unity (eq.~\ref{eq:psi}), and matrices
+$\mathcal{E}$ and $\mathcal{B}$ (eqs.~\ref{eq:matrix-E} and \ref{eq:matrix-B}).
+\end{itemize}
+
+
+
+
+
+
+
+
+
+
+
+
+
+%--------------------------------------------------------------------
+\subsection{First step: Injection} \label{chap:injection-step}
+%--------------------------------------------------------------------
+
+
+%------------------------------------------------
+\subsubsection{Injecting the Energy Density}
+%------------------------------------------------
+
+
+In the injection step, the radiation is gathered from radiating sources and injected into the gas.
+Radiation emitting sources are taken to be stars or entire stellar populations, which in SWIFT are
+represented by individual star particles. The equation to be solved in this step is
+
+\begin{equation*}
+ \DELDT{E_i} = \dot{E}_i^* \label{eq:solve-dEdt}
+\end{equation*}
+
+where $\dot{E}_i^*$ is the total energy density in the frequency group $i$ emitted by all stars $k$
+that are within compact support radius of the corresponding particle. If each star $k$ deposits
+some fraction $\xi_k$ of its respective total radiation energy density to be injected,
+then the total radiation energy density injected into a single gas particle $i$ is given by:
+
+\begin{align}
+ \dot{E_i}^* = \sum_k E_{i,k}^* \xi_k \label{eq:injection_onto_particle} \ .
+\end{align}
+
+The exact choice how the fraction $\xi_k$ is determined will be discussed further below.
+
+Eq.~\ref{eq:injection_onto_particle} is solved using a simple finite difference
+discretization and first order forward Euler integration:
+
+\begin{equation}
+ E_i(t + \Delta t) = E_i(t) + \Delta t \dot{E}_i^*
+\end{equation}
+
+for each particle.
+
+Star particles also have individual time step sizes, just like gas particles. Because both star
+and gas particles have individual time step sizes, the way the energy density $\dot{E}_i^*$ is
+deposited from stars onto particles is determined by the star particles' time steps. Each time a
+star particle is active, i.e. the simulation is at a step where the star particle finishes its time
+integration, the star particle ``pushes'' radiation onto both active and inactive gas particles.
+
+
+
+We now define the weights $\xi_s$ used in eq.~\ref{eq:injection_onto_particle} to distribute the
+total radiation energy ejected by a star $s$ over a time step $\Delta t$ onto the surrounding gas
+particles $p$.
+A natural way of depositing the energy density from the star on a neighboring gas particle is to
+make use of the already present partition of unity, $\psi_p(\x_p - \x_s)$ (see eqs. \ref{eq:psi} -
+\ref{eq:omega}), which is a fundamental building block for finite volume particle methods.
+% It guarantees to sum up to unity, and additionally is taking into account the particle
+% configuration due to its normalization.
+
+
+Then the total emitted energy of the star $s$, $E^*_i(\x_s)$, is fully and self-consistently
+distributed on the gas (here gas particles have the index $g$):
+
+\begin{align}
+ \sum_g \psi_g(\x_s)\ E_i^*(\x_s) = \sum_g \psi(\x_s - \x_g, h_s) \ E_i^*(\x_s) = E_i^*(\x_s)
+\end{align}
+
+For the actual distribution of the radiation from stars to gas particles, the sum over neighboring
+stars from the gas particle's point of view looks like:
+
+\begin{align}
+ E^*_i (\x = \x_{gas}) = \sum_{stars} E^*_i(\x_{star}) \ \psi(\x_{star} - \x_{gas}, h(\x_{star}))
+\end{align}
+
+suggesting
+
+\begin{align}
+ \xi_s = \psi(\x_{s} - \x_{g}, h(\x_{s})) \ .
+\end{align}
+
+
+A potential problem here is that unless the particle distribution is sufficiently symmetric, the
+resulting injected energy density (and flux) will not be isotropic, which is what stellar radiation
+is typically assumed to be. So instead, we ``correct'' the weight $\psi_p$ in an attempt to improve
+the isotropy of the deposited net energy and flux. Suppose we divide space into 8 octants centered
+on a star. Each octant $i$ is assigned a weight $w_i$ defined as
+
+\begin{align*}
+ w_i = \sum\limits_{\text{particles $j$ in $i$}} \xi (\x_j)
+\end{align*}
+
+We now apply the correction through a factor $\mu_i$ for each octant individually.
+Demanding that
+
+\begin{enumerate}
+\item The total weight must remain constant, i.e.
+$\sum_i \mu_i w_i = \sum_i w_i $
+
+\item The modified weight of each quadrant should be equal, such that at least along these eight
+directions, the injected energy and flux are isotropic:
+$\mu_i w_i = \mu_j w_j $ for all $i, j$
+
+\end{enumerate}
+
+gives us:
+
+\begin{equation}
+\mu_a = \frac{\sum_b w_b}{8 w_a}
+\end{equation}
+
+The correction term needs a small modification for cases where there is a octant which contains zero
+particles and hence zero weight: Let $q_{nz}$ be the number of octants with non-zero weight, i.e.
+with non-zero particles in them. Then
+
+\begin{equation}
+ \mu_a = \frac{\sum_b w_b}{q_{nz} w_a} \label{eq:isotropy-correction-with-zero}
+\end{equation}
+
+
+this ultimately gives us
+
+\begin{align}
+ \xi_s = \mu_a \ \psi(\x_{s} - \x_{g}, h(\x_{s}))
+\end{align}
+
+where $a$ is the octant in which the gas particle resides.
+
+
+In principle, we could also inject photon fluxes directly from stellar sources, just as we do with
+photon energy density. Test however have shown that best results are achieved by not injecting any
+net flux onto gas particles, so we don't do that in GEAR-RT.
+
+
+
+
+
+
+
+
+
+
+
+
+%-----------------------------------------------------------------------
+\subsection{Second Step: Transport}\label{chap:transport-step}
+%-----------------------------------------------------------------------
+
+%-----------------------------------------------------------------------
+\subsubsection{Solving the Transport Equations}
+%-----------------------------------------------------------------------
+
+In this step, we solve
+\begin{align}
+ \DELDT{E_i} + \nabla \cdot \Fbf_i &= 0 \label{eq:homogenized-dEdt}\\
+ \DELDT{\Fbf_i} + c^2 \nabla \cdot \mathds{P}_i &= 0 \label{eq:homogenized-dFdt}
+\end{align}
+
+They are both in the form of a hyperbolic conservation law, so we can use the finite volume
+particle method given in eq.~\ref{eq:meshless-explicit}.
+
+We have the state vector $\U$ and flux tensor $\F$:
+
+\begin{align}
+ \U =
+ \begin{pmatrix}
+ E_i \\
+ \Fbf_i
+ \end{pmatrix}
+ \quad , \quad
+ \F =
+ \begin{pmatrix}
+ \Fbf_i \\
+ c^2 \mathds{P}_i
+ \end{pmatrix}
+\end{align}
+
+The total sum $\sum_l \F_{kl} A_{kl}$ for each particle $k$ is collected during a neighbor
+interaction loop. Once the sum is accumulated, the final explicit update of the state is given by:
+
+\begin{align}
+\U_k(t + \Delta t_k) =
+ \U_k (t) - \frac{1}{V_k} \sum\limits_{l} \min\{\Delta t_k, \Delta t_l\}
+ \F_{\alpha,kl} \mathcal{A}_{kl}^\alpha
+\label{eq:transport-update-explicit}
+\end{align}
+
+$V_k$ are the associated particle volumes given in eq.~\ref{eq:volume}. \Aij are the effective
+surface terms given in eq.~\ref{eq:HopkinsAij}.
+
+
+
+
+
+It remains to determine the inter-particle fluxes $\F_{kl}$. We also follow the scheme outlined by
+the finite volume particle method. As it is done for the hydrodynamics, we approximate the problem
+by defining an ``interface'' at the position
+
+\begin{equation}
+ \mathbf{x}_{kl} = \mathbf{x}_k + \frac{h_k}{h_k + h_l} ( \mathbf{x}_l - \mathbf{x}_k )
+\end{equation}
+
+and extrapolate the value of the conserved variables at the position:
+
+\begin{equation}
+ \U_{k, l}(\mathbf{x} =
+ \mathbf{x}_{kl}) \approx \U_{k, l} + \DELDX{ \U_{k, l}}\ (\mathbf{x}_{kl}
+- \mathbf{x}_{k,l}) \label{eq:gradient-extrapolation}
+\end{equation}
+
+
+such that the fluxes $\F_{kl}$ in the update formula~\ref{eq:transport-update-explicit} between
+particles is estimated as the solution of the Riemann problem with ``left'' state $\U_k(\x_{kl})$
+and ``right'' state $\U_l(\x_{kl})$
+
+\begin{align}
+ \F_{kl}
+ = RP
+ \left(\U_k + (\x_{k} - \x_{kl})_\alpha \DELDXALPHA{\U_k}, \
+ \U_l + (\x_{l} - \x_{kl})_\alpha \DELDXALPHA{\U_l} \right) \label{eq:rt-riemann-problem}
+\end{align}
+
+
+GEAR-RT offers two approximate Riemann solvers to find solutions to this Riemann problem. They
+are described in Section~\ref{chap:riemann-rt}.
+
+The gradients $\DELDX{\U_{k,l}}$ are estimated again using the least-squares second order accurate
+gradient expression given in eqn. \ref{eq:gradient}, for which a particle-particle interaction loop
+is required, which furthermore needs to be performed before the transport interaction loop, during
+which the sum of the exchanged (time integrated) fluxes of eq.~\ref{eq:transport-update-explicit}
+are accumulated. Since the RT is solved after the hydrodynamics and the particle positions haven't
+changed since the last hydrodynamics drift, the particles' smoothing lengths are still accurate, and
+we don't require an additional density loop for the radiative transfer. For that reason the gradient
+particle loop constitutes the first particle interaction loop in the RT scheme.
+
+
+
+
+%---------------------------------------------------
+\subsubsection{Flux Limiters}
+%---------------------------------------------------
+
+The gradient extrapolation (eq.~\ref{eq:gradient-extrapolation}) is equivalent from going from a
+piece-wise constant reconstruction of the radiation field to a piece-wise linear one, which makes
+the method second order accurate in space. For each quantity $Q_{k,l} \in \U = (E_{k,l},
+\Fbf_{k,l})^T$, we extrapolate the value at the interface $\x_{kl}$ using the estimated gradient
+$dQ_{k,l}$. To prevent oscillations and instabilities from occurring, we need to employ a flux
+limiter at this point, which effectively reduces the gradient slopes $dQ_k$ and $dQ_l$ whenever
+necessary.
+
+
+The two limiters that worked well in tests are the \emph{minmod limiter}
+
+\begin{align}
+% \alpha_{\text{minmod}}(dQ_k, dQ_l) =
+% \begin{cases}
+% \text{if } dQ_k \times dQ_l > 0:
+% \begin{cases}
+% \text{if } |dQ_k| < |dQ_l| :\quad &\alpha = dQ_k / dQ_l \\
+% \text{else: } &\alpha = dQ_l / dQ_k\\
+% \end{cases}\\
+% \text{else: } \alpha = 0
+% \end{cases}
+ dQ_{\text{minmod}}(dQ_k, dQ_l) =
+ \begin{cases}
+ dQ_k & \text{ if } |dQ_k| < |dQ_l| \text{ and } dQ_k dQ_l > 0 \\
+ dQ_l & \text{ if } |dQ_k| > |dQ_l| \text{ and } dQ_k dQ_l > 0 \\
+ 0 & \text{ if } dQ_k dQ_l \leq 0
+ \end{cases} \label{eq:rt-minmod}
+\end{align}
+
+and the \emph{monotonized central difference (MC) limiter}:
+\begin{align}
+ r &= dQ_k / dQ_l \nonumber \\
+ \alpha_{\text{MC}}(r) &= \max \left\{ 0, \min\left[\frac{1 + r}{2}, 2, 2r \right] \right\}
+ \label{eq:rt-MC}\\
+ dQ_{k, limited} &= \alpha_{MC}(r)\ dQ_k \nonumber\\
+ dQ_{l, limited} &= \alpha_{MC}(r)\ dQ_l \nonumber
+\end{align}
+
+
+The minmod limiter appears to be the most robust one with regards to ensuring stability for
+radiative transfer, and for this reason is the recommended choice.
+
+
+Other possible choices are the \emph{superbee} limiter:
+
+\begin{align}
+ r &= dQ_k / dQ_l \nonumber \\
+ \alpha_{\text{superbee}}(r) &= \max \left\{ 0, \min (1, 2r), \min(2, r) \right\}
+ \label{eq:rt-superbee}\\
+ dQ_{k,l,\ limited} &= \alpha_{superbee}(r)\ dQ_{k,l} \nonumber
+\end{align}
+
+
+And the \emph{van Leer} limiter:
+
+\begin{align}
+ r &= dQ_k / dQ_l \nonumber\\
+ \alpha_{\text{vanLeer}}(r) &= \frac{r + |r|}{1 + |r|}
+ \label{eq:rt-vanLeer}\\
+ dQ_{k,l,\ limited} &= \alpha_{vanLeer}(r)\ dQ_{k,l} \nonumber
+\end{align}
+
+
+
+
+
+
+
+
+%-------------------------------------------------------------------------------------------
+\subsubsection{Riemann Solvers for the Moments of the RT equation}\label{chap:riemann-rt}
+%-------------------------------------------------------------------------------------------
+
+It remains to find a solution for the Riemann problem~\ref{eq:rt-riemann-problem} which gives us
+the inter-particle (``inter-cell'') fluxes $\F_{kl}$. To find $\F_{kl}$, we use the extrapolated and
+flux limited states $E_{k,l}$ and $\Fbf_{k,l}$ at the interface position $\x_{kl}$ to compute the
+states and fluxes of the conservation law, $\U_{k,l}(\x_{kl})$ and $\F_{k,l}(\x_{kl})$.
+
+The components of the vector valued photon flux density $\Fbf_{k,l}$ are projected along the unit
+vector pointing from the particle towards the surface, allowing us to treat each component
+individually as a one dimensional problem. We then adapt the convention that particle $k$ is the
+left state of the Riemann problem, $\U_L$ and $\F_R$, while particle $l$ is the right state $\U_R$
+and $\F_R$.
+
+The \emph{Global Lax Friedrich (GLF)} Riemann solver \citep{ramses-rt13} then gives the
+following solution for the flux $\F$ at an interface which separates a left state $\U_L$
+and right state $\U_R$:
+\begin{equation}
+ \F_{kl}(\U_L, \U_R) =
+ \frac{\F_{L} + \F_{R}}{2} -
+ \frac{c}{2} \left(\U_R - \U_L \right) \label{eq:riemann-GLF}
+\end{equation}
+
+which gives us an approximate solution for the flux $\F_{kl}$ we are looking for.
+
+
+
+An alternative is the \emph{Harten - Lax - van Leer (HLL)} Riemann solver
+(\cite{gonzalezHERACLESThreedimensionalRadiation2007}), given by:
+
+\begin{align}
+ \F(\U_L, \U_R) &=
+ \frac{ \lambda^{+} \F_{L} - \lambda^{-} \F_{R} + \lambda^{+}
+\lambda^{-} (\U_R - \U_L)}{ \lambda^{+} - \lambda^{-} }
+\label{eq:riemann-HLL} \\
+ \lambda^+ &= \max(0, \lambda_{max})\\
+ \lambda^- &= \min(0, \lambda_{min})
+\end{align}
+
+
+Here, $\lambda^{max}$ and $\lambda^{min}$ are the minimum and the maximum of the Eigenvalues of the
+Jacobian matrix $\frac{\del \F(\U)}{\del \U}$. It turns out that the Eigenvalues can be described
+using only two parameters, the angle between the flux and the interaction interface, $\theta$, and
+the reduced flux $\mathbf{f} = \Fbf / (cE)$. Rather then continuously computing them on-the-fly for
+every interaction, they are tabulated\footnote{
+A program to produce the tables of the Eigenvalues depending on $\theta$ and $\mathbf{f}$ is
+provided on \url{https://github.com/SWIFTSIM/swiftsim-rt-tools}.
+}
+and interpolated.
+
+The interest in the HLL solver is mainly because while it is more expensive than the GLF solver, it
+was shown \citep{ramses-rt13, gonzalezHERACLESThreedimensionalRadiation2007} to be less diffusive
+and is thus better suited to form shadows more correctly than the GLF solver, which is a known
+weakness of the moment based radiative transfer.
+
+
+
+
+
+
+
+
+
+
+
+%---------------------------------------------
+\subsection{Third step: Thermochemistry}
+%---------------------------------------------
+
+
+In this final step, we solve for the interaction between photons and the gas. The equations to be
+solved are:
+
+\begin{align*}
+ \DELDT{E_i} &=
+ - \sum\limits_{j}^{\absorbers} n_j \sigma_{i j} c E_i + \dot E_i ^{rec}
+\\
+ \DELDT{\Fbf_i} &=
+ - \sum\limits_{j}^{\absorbers} n_j \sigma_{i j} c \Fbf_i
+\end{align*}
+
+The source term for recombination radiation, $\dot{E_i}^{rec}$, was added in for completeness.
+GEAR-RT currently (2023) neglects the emission of recombination.
+
+The actual thermochemistry is solved using the GRACKLE library\footnote{
+GRACKLE is available on \url{https://github.com/grackle-project/grackle}.\\
+Note that a frozen version of GRACKLE, tested and verified to work with GEAR-RT and other SWIFT
+projects, is being maintained on \url{https://github.com/mladenivkovic/grackle-swift}. It is
+strongly suggested to use that version of GRACKLE.
+}
+\citep{smithGrackleChemistryCooling2017}. This involves the evolution of the individual species
+number densities, as well as the heating of the gas. Currently (2023) GEAR-RT supports the ``6
+species network'', composed of H$^0$, H$^+$, He$^0$, He$^+$, He$^{++}$, and $e^{-}$, provided by
+GRACKLE.
+
+GRACKLE requires us to provide it with the (total) photo-heating rate $\mathcal{H} =
+\sum_{i,j} \mathcal{H}_{ij}$ and the individual photo-ionization rates $\Gamma_{j} = \sum_i
+\Gamma_{i,j}$ for each species $j$ and photon group $i$. The expressions for the heating and
+photo-ionization rates for each photon group are given in eqns.~\ref{eq:photoheating-group} and
+\ref{eq:photoionization-group}. GRACKLE then solves the network of equations of thermochemistry for the
+6 species and evolves the internal energy of the gas over some time step
+$\Delta t$.
+
+While GRACKLE evolves the state of the gas, we need to take care of the removal of the absorbed
+radiation energy over the thermochemistry time step manually. We solve this again using a simple
+first order forward Euler integration. A minor complication is that the number density of the
+ionizing species $n_j$ is not constant over the time step $\Delta t$, as we are currently in the
+process of ionizing the gas. This may lead to scenarios where the absorption rate is overestimated
+by not accounting for the ionization during the thermochemistry time step. To account for this
+effect, we take the average absorption rate over the time step instead:
+
+
+\begin{align*}
+ \frac{E_i (t + \Delta t) - E_i}{\Delta t}
+ = -\overline{\epsilon}_i \ N_i(t) c \ \sum_j \sigma_{i j} ^ N \
+ \frac{n_j(t + \Delta t) + n_j(t)}{2}
+\end{align*}
+
+which we can do since GRACKLE provides us with the updated value $n_j(t + \Delta t)$.
+The final equation used to remove absorbed photons from the radiation field is:
+
+\begin{align}
+E_i (t + \Delta t)
+ = E_i(t) \left(1 - c \ \sum_j \sigma_{i j} ^ N \
+ \frac{n_j(t + \Delta t) + n_j(t)}{2} \ \Delta t \right)
+\end{align}
+
+The same is applied to the photon fluxes:
+
+\begin{align}
+\Fbf_i (t + \Delta t)
+ &= \Fbf_i(t) \left(1 - c \ \sum_j \sigma_{i j} ^ N \
+ \frac{n_j(t + \Delta t) + n_j(t)}{2} \ \Delta t \right)
+\end{align}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+%---------------------------------------------
+\subsection{Additional Topics}
+%---------------------------------------------
+
+
+%-----------------------------------------------------
+\subsubsection{Computing the Time Step Size}
+%-----------------------------------------------------
+
+The time step size of the particles is computed in the same manner as we do it for hydrodynamics.
+A ``cell size'' $\Delta x$ of a particle is estimated using
+
+\begin{align}
+ \Delta x \approx \left(\frac{V_i}{V_{S,\nu}} \right)^{1/\nu}
+\end{align}
+
+where $V_{S,\nu}$ is the volume of a $\nu$-dimensional unit sphere, i.e. $4 \pi / 3$ in 3D. The
+signal velocity in case of radiation is simply the speed of light $c$, which leads to the time step
+estimate
+
+\begin{align}
+ \Delta t = C_{CFL} \frac{\Delta x}{c} \label{eq:rt-cfl}
+\end{align}
+
+
+
+%---------------------------------------------
+\subsubsection{Dealing With Particle Drifts}\label{chap:rt-drift}
+%---------------------------------------------
+
+While we assume particles are static w.r.t. the simulated volume in the context of RT, they are
+drifted for the purposes of hydrodynamics. As a consequence, we need to make corrections to the
+radiation fields when a particle is drifted for hydrodynamics purposes. We do that by simply
+extrapolating the particle value at the new point given its current state and the gradients $\nabla
+\U$ we obtained using the general gradient expression given in eq~\ref{eq:gradient}. Explicitly, for
+a particle $i$ and for each component of the radiation state vector $\U^k$, we first obtain the
+drift distance over the time $\Delta t_{i, drift}$ as
+
+\begin{align}
+\Delta \x_{i,\text{drift}} &= \V_{i} \Delta t_{i,\text{drift}}
+\end{align}
+
+where $\V_i$ is the particle velocity. We then obtain the value corrected at the drifted
+location
+
+\begin{align}
+\U_{i}^k (\x + \Delta \x_{\text{drift}}) \approx \U_{i}^k (\x) + \nabla \U_{i}^k (\x) \cdot \Delta
+x_{i,\nu,\text{drift}}
+\end{align}
+
+
+
+
+
+
+
+
+
+%---------------------------------------------
+\subsubsection{Cosmology}
+%---------------------------------------------
+
+TODO
+
+
+
+
+
+
+
+
+%---------------------------------------------
+\subsubsection{Reducing The Speed of Light}
+%---------------------------------------------
+
+
+In an attempt to reduce the computational expense associated with these tiny time step sizes,
+\cite{ramses-rt13} employ the Reduced Speed of Light Approximation (RSLA), which was introduced by
+\cite{gnedinMultidimensionalCosmologicalRadiative2001}. I also adapt the RSLA with GEAR-RT. The
+approach consists of globally reducing the speed of light everywhere by some user-set factor $f_c$,
+i.e. replace the speed of light $c$ in all equations involved with radiative transfer by
+
+\begin{align}
+ \tilde{c} = f_c c
+\end{align}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+%---------------------------------------------
+\subsubsection{Ion Mass Fluxes}
+%---------------------------------------------
+
+Hydrodynamics with Finite Volume Particle Methods exchange mass fluxes $\Delta m$ between particles.
+This means that we need to pay attention to the individual ionizing and ionized species' mass
+fractions being exchanged. The mass of each individual species being exchanged needs to be traced
+individually, and depending on the direction: If a particle loses mass, then it loses mass and
+species mass fractions according to its own mass fractions. If it gains mass during an exchange,
+then it gains species mass fractions according to the interaction partner's mass fractions.
+
+The usual convention for a mass flux between particles ``$i$'' and ``$j$'' is to treat $i$ as the
+``left'' particle and $j$ as the ``right'' particle in an analogy to cells. Let $X_k$ denote the
+mass fraction of each species $k$. During each \emph{hydrodynamics} flux exchange, for each species
+$k$ the total mass flux $\Delta x_{i,k}$ is being accumulated:
+
+\begin{align}
+\Delta x_{i,k} &= \sum_j \Delta m_{ij} X_{ij,k}\\
+\text{where }\quad X_{ij,k} &=
+ \begin{cases}
+ X_{i,k} \quad \text{ if } \quad \frac{\Delta m_{ij}}{\Delta t} > 0 \\
+ X_{j,k} \quad \text{ if } \quad \frac{\Delta m_{ij}}{\Delta t} < 0
+ \end{cases}
+\end{align}
+
+The masses which particles carry are updated during the kick operation, during which the individual
+species' masses $x_{i,k}$ are also evolved. For a particle $i$ with mass $m_i$, the masses of
+individual species at the beginning of a step are given by:
+
+\begin{align}
+ x_{i,k}^{t} &= m_i X_{i,k}^{t}
+\end{align}
+
+During the kick operation, they get updated as follows:
+
+\begin{align}
+ x_{i,k}^{t + \Delta t} &= x_{i,k}^{t} + \Delta x_{i,k}
+\end{align}
+
+Finally, after the update, the new mass fractions of the species can be obtained using
+
+\begin{align}
+X_{i,k}^{t + \Delta t} &= \frac{x_{i,k}^{t + \Delta t}}{\sum_k x_{i,k}^{t + \Delta t}}
+\end{align}
+
+
+
diff --git a/theory/RadiativeTransfer/GEARRT/GEARRT.tex b/theory/RadiativeTransfer/GEARRT/GEARRT.tex
new file mode 100644
index 0000000000000000000000000000000000000000..e90a2aad2a0524624763079e47f326395a4c66fe
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/GEARRT.tex
@@ -0,0 +1,57 @@
+\input{./header}
+
+
+
+%------------------------------------------
+%: Metadata
+%------------------------------------------
+
+\title{GEAR-RT}
+\subtitle{Theory Documentation}
+\author{Mladen Ivkovic}
+\date{September 2023}
+
+%------------------------------------------
+
+
+
+
+
+
+\begin{document}
+
+
+%Titlepage
+\maketitle
+\tableofcontents
+\clearpage
+
+
+
+
+
+\include{1-meshless}
+\include{2-rt-equations}
+\include{3-rt-numerical-method}
+
+
+
+
+
+
+
+
+\clearpage
+\appendix{}
+\include{appendix}
+
+
+%------------
+%Bibliography
+%------------
+\clearpage
+\bibliography{references}
+
+
+
+\end{document}
diff --git a/theory/RadiativeTransfer/GEARRT/appendix.tex b/theory/RadiativeTransfer/GEARRT/appendix.tex
new file mode 100644
index 0000000000000000000000000000000000000000..b739b91ac701c983549c20e6a54548671d4fb381
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/appendix.tex
@@ -0,0 +1,288 @@
+%------------------------------------------------
+\section{RT Related Quantities And Their Units}\label{app:variables}
+%------------------------------------------------
+
+It could be helpful to have the commonly used quantities of RT written down
+somewhere along with their units. So here you go.
+
+
+\include*{./tables/rt_variables}
+
+
+
+
+
+
+
+\clearpage
+%--------------------------------------------
+\section{Function Calls of an RT step}
+%--------------------------------------------
+
+
+
+\begin{landscape}
+{\footnotesize
+
+\begin{tabular}[l]{%
+ >{\raggedright\arraybackslash}p{2.6cm}%
+ >{\raggedright\arraybackslash}p{2.8cm}%
+ >{\raggedright\arraybackslash}p{7cm}%
+ >{\raggedright\arraybackslash}p{7cm}%
+}
+\textbf{task} & \textbf{task type} & \textbf{task purpose} & \textbf{function calls} \\[.5em]
+\hline
+\hline
+Injection Prep &
+ Interacting star and gas particles &
+ Gather gas particle neighbour data in preparation for the injection &
+ \texttt{runner\_iact\_rt\_inject\_prep} in \verb|src/rt/method/rt_iact.h| \\
+\hline
+Star Emission Rates &
+ Work on individual star particles &
+ Prepare everything necessary that needs to be done for radiation sources before the radiation sources and the gas interact. &
+ \texttt{rt\_compute\_stellar\_emission\_rate} in \verb|src/rt/method/rt.h| \\
+\hline
+\hline
+RT in &
+ Implicit&
+ Collect dependencies &
+ - \\
+\hline
+Injection &
+ Interacting star and gas particles &
+ Distribute the radiation from star particles onto gas particles &
+ \verb|runner_iact_rt_inject| in \verb|src/rt/method/rt_iact.h| \\
+\hline
+Ghost1 &
+ Work on individual gas particles &
+ Finish up any work that needs to be done for gas particles before the next gas $\leftrightarrow$ gas interaction begins &
+ \texttt{rt\_finalise\_injection} in \verb|src/rt/method/rt.h|\\
+\hline
+Gradient &
+ Interacting gas and gas particles &
+ Compute necessary gradients of the radiation quantities &
+ \verb|rt_gradients_collect| in \verb|src/rt/method/rt_gradients.h| and
+ \verb|rt_gradients_nonsym_collect| in \verb|src/rt/method/rt_gradients.h|\\
+\hline
+Ghost2 &
+ Work on individual gas particles &
+ Finish up any work that needs to be done for gas particles before the next gas $\leftrightarrow$ gas interaction begins &
+ \texttt{rt\_end\_gradient} in \verb|src/rt/method/rt.h|\\
+\hline
+Transport &
+ Interacting gas and gas particles &
+ Compute/exchange fluxes of homogenized equation of radiative transfer. &
+ \verb|runner_iact_rt_transport| in \verb|src/rt/method/rt_iact.h| and
+ \verb|runner_iact_nonsym_rt_transport| in \verb|src/rt/method/rt_iact.h|\\
+\hline
+Transport out &
+ Implicit&
+ Collect dependencies &
+ - \\
+\hline
+Thermochemistry &
+ Work on individual gas particles &
+ Finish up any work that needs to be done for gas particles before the thermochemistry part of the computation can begin. Then do the thermochemistry computation. &
+ \texttt{rt\_finalise\_transport} in \verb|src/rt/method/rt.h|,
+ \verb|rt_do_thermochemistry| in \verb|src/rt/method/rt_thermochemistry.h|\\
+\hline
+RT out &
+ Implicit&
+ Collect dependencies &
+ - \\
+\hline
+\end{tabular}
+
+
+}
+\end{landscape}
+
+
+
+
+
+
+
+
+
+\clearpage
+%-------------------------------------------------------------------------
+\section{Creating Collisional Ionization Equilibrium Initial Conditions}
+%-------------------------------------------------------------------------
+
+
+On startup, GEARRT offers the possibility to generate the ionization mass fractions of the gas
+particles assuming the gas is in collisional ionization equilibrium, composed of Hydrogen and
+Helium, and that there is no radiation present. In order to determine the ionization mass fractions
+of all species (H$^0$, H$^+$, He$^0$, He$^+$, He$^{++}$) for a given specific internal energy $u$,
+an iterative procedure needs to be applied because the gas variables are interconnected:
+
+
+
+Commonly the average particle mass $\overline{m}$ of a gas is expressed using the (unitless) mean
+molecular weight $\mu$ :
+
+\begin{align}
+ \overline{m} = \mu m_u
+\end{align}
+
+where $m_u$ is the atomic mass unit and $\mu$ is given by
+
+\begin{align}
+ \frac{1}{\mu} = \sum_j \frac{X_j}{A_j} (1 + E_j)
+\end{align}
+
+
+Specifically, ionization changes the mass fractions $X_j$ of the species, and therefore also the
+mean molecular weight $\mu$. In turn, the mean molecular weight determines the gas temperature at a
+given specific internal energy:
+
+\begin{align}
+ T = u (\gamma - 1) \mu \frac{m_u}{k}
+\end{align}`
+
+Lastly, the gas temperature determines the collisional ionization and recombination rates, which
+need to be balanced out by the correct number density of the individual species in order to be in
+ionization equilibrium, i.e. at a state where the ionization and recombination rates exactly cancel
+each other out.
+
+We take the ionization and recombination rates from \citet{katzCosmologicalSimulationsTreeSPH1996},
+which are given in Table \ref{tab:coll-ion-rates-katz}. For a gas with density $\rho$, Hydrogen
+mass fraction $X_H$ and Helium mass fraction $X_{He} = 1 - X_H$, the total number densities of all
+hydrogen and helium species are
+
+\begin{align}
+ n_H &= X_H \frac{\rho}{m_u} \\
+ n_{He} &= X_{He} \frac{\rho}{4 m_u}
+\end{align}
+
+and in equilibrium, the number densities of the individual species are given by
+
+\begin{align}
+ n_{H^0} &=
+ n_H \frac{A_{H^+}}{A_{H^+} + \Gamma_{H^0}} \\
+ n_{H^+} &=
+ n_H - n_{H^0} \\
+ n_{He^+} &=
+ n_{He} \frac{1}{1 + (A_{He^+} + A_d) / \Gamma_{He^0} + \Gamma_{He^+} / A_{He^{++}}} \\
+ n_{He^0} &=
+ n_{He^+} \frac{A_{He^+} + A_d}{\Gamma_{He^0}} \\
+ n_{He^{++}} &=
+ n_{He^+} \frac{\Gamma_{He^+}}{A_{He^+}}
+\end{align}
+
+
+To summarize, the tricky bit here is that the number densities determine the mean molecular weight,
+the mean molecular weight determines the temperature of the gas for a given density and specific
+internal energy, while the temperature determines the number densities of the species.
+
+To find the correct mass fractions, the iterative Newton-Raphson root finding method is used.
+Specifically, using some initial guesses for temperature and mean molecular weights, $T_{guess}$
+and $\mu_{guess}$, in each iteration step we determine the resulting specific internal energy
+$u_{guess} = k T_{guess} / (\gamma - 1) / (\mu_{guess} m_u)$. The function whose root we're looking
+for is
+
+\begin{align}
+ f(T) = u - u_{guess}(T)
+\end{align}
+
+with the derivative
+
+\begin{align}
+ \frac{\del f}{\del T} = - \frac{\del u}{\del T} (T = T_{guess}) = \frac{k}{(\gamma -
+1) / (\mu_{guess} m_u )}
+\end{align}
+
+where $u$ is the specific internal energy of the gas, which is fixed and provided by the initial
+conditions. We now look for the $T$ at which $f(T) = 0$. The Newton-Raphson method prescribes to
+find the $n+1$th $T_{guess}$ using
+
+\begin{align}
+ T_{n+1} = T_n + \frac{f(T_n)}{\frac{\del f}{\del T}(T_n)}
+\end{align}
+
+During each iteration, the new mass fractions and the resulting mean molecular weight given the
+latest guess for the temperature are computed. At the start, the first guess for the temperature
+$T_{guess}$ is computed assuming a fully neutral gas. Should that gas temperature be above $10^5$ K,
+the first guess is changed to a fully ionized gas. The iteration is concluded once $f(T) \leq
+\epsilon = 10^{-4}$.
+
+
+
+\include*{tables/thermochemistry_rates}
+
+
+
+
+
+
+
+
+
+\clearpage
+%-------------------------------------------------------------------------------
+\section{Converting Photon Number Emission Rates to Photon Energy Emission Rates}
+%-------------------------------------------------------------------------------
+
+Many other, in particular older, codes and papers use photon \emph{number} injection rates $\dot{N}_{\gamma}$ for their emission rates rather than the \emph{energy} injection rates $\dot{E}_\gamma$ or equivalently luminosities $L$.
+
+To convert between these two quantities, we need to assume that the emission follows some spectrum $J(\nu)$.
+
+In the case of a single photon group, the conversion is quite simple: We first need to compute the average photon energy $\overline{E}_\gamma$:
+
+\begin{align*}
+ \overline{E}_\gamma = \frac{\int J(\nu) \ \de \nu }{\int J(\nu) / (h \nu) \ \de \nu}
+\end{align*}
+
+then the emitted luminosity (energy per unit time) is
+
+\begin{align}
+ L = \overline{E}_\gamma \ \dot{N}_{\gamma}
+\end{align}
+
+Note that in many cases, the given emission photon number rate is the number rate of \emph{ionizing} photons. For us, this means that we need to start the integrals at the lowest ionizing frequency $\nu_{\text{ion, min}}$ in order to have the correct translation to the luminosity of the \emph{ionizing} energy:
+
+
+\begin{align}
+ \overline{E}_\gamma = \frac{\int\limits_{\nu_{\text{ion, min}}}^\infty J(\nu) \ \de \nu }{\int\limits_{\nu_{\text{ion, min}}}^\infty J(\nu) / (h \nu) \ \de \nu}
+\end{align}
+
+
+In the case of several photon groups being used, the conversion requires a little adaptation in order to preserve the correct number of photons emitted. For each photon group $i$, the average photon energy is given by
+
+\begin{align}
+ \overline{E}_i = \frac{\int\limits_{\nu_{i \text{, min}}}^{\nu_{i \text{, max}}} J(\nu) \ \de \nu }{\int\limits_{\nu_{i \text{, min}}}^{\nu_{i \text{, max}}} J(\nu) / (h \nu) \ \de \nu}
+\end{align}
+
+
+Secondly, we need to compute the fraction $f_i$ of ionizing photons in each bin, which is given by
+
+\begin{align}
+ f_i = \frac{\int\limits_{\nu_{i \text{, min}}}^{\nu_{i \text{, max}}} J(\nu) / (h \nu) \ \de \nu }{\int\limits_{\nu_{min}}^{\infty} J(\nu) / (h \nu) \ \de \nu}
+\end{align}
+
+Then the number of emitted photons in each bin is given by
+
+\begin{align}
+\dot{N}_i = f_i\ \dot{N}_\gamma
+\end{align}
+
+And the luminosities are given by
+
+\begin{align}
+ L_i &= \overline{E}_i \ \dot{N}_i \\
+ &= \frac{\int\limits_{\nu_{i \text{, min}}}^{\nu_{i \text{, max}}} J(\nu) \ \de \nu }{\int\limits_{\nu_{min}}^{\infty} J(\nu) / (h \nu) \ \de \nu} \ \dot{N}_\gamma
+\end{align}
+
+
+
+Python scripts to compute and convert photon number rates into luminosities are provided in
+\url{https://github.com/SWIFTSIM/swiftsim-rt-tools}.
+
+
+
+
+
+
+
diff --git a/theory/RadiativeTransfer/GEARRT/header.tex b/theory/RadiativeTransfer/GEARRT/header.tex
new file mode 100644
index 0000000000000000000000000000000000000000..4d70ae5486be35478305da801139eca81a3b59b1
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/header.tex
@@ -0,0 +1,204 @@
+\documentclass[12pt, a4paper, english, singlespacing, parskip]{scrartcl}
+
+%\documentclass[
+%11pt, % The default document font size, options: 10pt, 11pt, 12pt
+%oneside, % Two side (alternating margins) for binding by default, uncomment to switch to one side
+%chapterinoneline, % Have the chapter title next to the number in one single line
+%english, % ngerman for German
+%singlespacing, % Single line spacing, alternatives: onehalfspacing or doublespacing
+%draft, % Uncomment to enable draft mode (no pictures, no links, overfull hboxes indicated)
+%nolistspacing, % If the document is onehalfspacing or doublespacing, uncomment this to set spacing in lists to single
+%liststotoc, % Uncomment to add the list of figures/tables/etc to the table of contents
+%toctotoc, % Uncomment to add the main table of contents to the table of contents
+%parskip, % Uncomment to add space between paragraphs
+%nohyperref, % Uncomment to not load the hyperref package
+%headsepline, % Uncomment to get a line under the header
+%]{scrartcl or scrreprt or scrbook} % The class file specifying the document structure
+
+\usepackage{lmodern} % Diese beiden packages sorgen für echte
+\usepackage[T1]{fontenc} % Umlaute.
+
+\usepackage{amssymb, amsmath, color, graphicx, float, setspace, tipa}
+\usepackage[utf8]{inputenc}
+\usepackage[english]{babel}
+\usepackage[pdfpagelabels,
+ pdfstartview = FitH,
+ bookmarksopen = true,
+ bookmarksnumbered = true,
+ linkcolor = black,
+ plainpages = false,
+ hypertexnames = false,
+ citecolor = black,
+ breaklinks]{hyperref}
+\usepackage{url}
+\usepackage{array}
+
+\allowdisplaybreaks % allows page breaks in align/equation environment
+
+\usepackage{authblk} % titlepage stuff
+\usepackage[titletoc, title]{appendix}
+
+\usepackage{newclude} % use \include*{file} instead \include{} to omit pagebreak after include
+\usepackage{lscape}
+
+
+
+
+%===================
+% BIBLIOGRAPHY
+%===================
+
+\usepackage[]{natbib}
+\bibliographystyle{apalike}
+
+%DONT FORGET TO COMPILE THE BIBLIOGRAPHY WITH BIBTEX WHEN CHANGES ARE MADE.
+
+
+
+%--------------------------------------------
+% OPTIONAL
+%--------------------------------------------
+
+
+
+%% Change font
+%\newcommand{\changefont}[3]{
+%\fontfamily{#1} \fontseries{#2} \fontshape{#3} \selectfont}
+%\changefont{ppl}{m}{n} nach \begin{document} einsetzen
+
+% Fig. instead of Figure, Tab. instead of Table
+%\usepackage[footnotesize]{caption2}
+%\addto\captionsenglish{\renewcommand{\figurename}{Fig.}}
+%\addto\captionsngerman{\renewcommand{\figurename}{Fig.}}
+%\renewcommand{\tablename}{Tab.}%
+
+%\pagestyle{headings} % Write headings on each page
+
+%\usepackage{chngcntr} \counterwithout{figure}{section} % Integer only figure numbers, ignoring chapter numbers
+
+
+
+
+
+
+%-----------------------------------------------
+% FORMAT TITLE
+%-----------------------------------------------
+
+
+% Set fonts of document parts
+\setkomafont{title}{\rmfamily\bfseries\boldmath}
+\addtokomafont{section}{\rmfamily\bfseries\boldmath}
+\addtokomafont{subsection}{\rmfamily\bfseries\boldmath}
+\addtokomafont{subsubsection}{\rmfamily\bfseries\boldmath}
+\addtokomafont{disposition}{\rmfamily} % table of contents and stuff
+\setkomafont{descriptionlabel}{\rmfamily\bfseries\boldmath}
+
+
+\usepackage{dsfont} % for Pressure tensor P with \mathds{}
+
+
+
+
+%-----------------------------------------------
+% Document specific definitions
+%-----------------------------------------------
+
+\newcommand{\Aij}{$\mathcal{A}_{ij}$} % A_ij
+\newcommand{\Aijm}{\ensuremath{\mathcal{A}_{ij}}} % A_ij math
+\newcommand{\U}{\ensuremath{\mathcal{U}}} % State vector
+\newcommand{\W}{\ensuremath{\mathcal{W}}} % State vector
+\newcommand{\F}{\ensuremath{\mathcal{F}}} % Flux tensor
+\newcommand{\Ubf}{\ensuremath{\mathbf{U}}} % State vector
+\newcommand{\Fbf}{\ensuremath{\mathbf{F}}} % Flux tensor
+\newcommand{\psitilde}{\ensuremath{\tilde{\psi}}} % psi tilde
+\newcommand{\half}{1/2} % 1/2
+
+
+\newcommand{\absorbers}{\text{HI, HeI, HeII}}
+\newcommand{\Ndot}{\dot{N}}
+
+
+%-----------------------------------------------
+% General Lazyness
+%-----------------------------------------------
+
+\newcommand{\corresponds}{\mathrel{\widehat{=}}} % equals with hat
+
+\newcommand {\arctanh}{\mathrm{arctanh}} % Atanh
+\newcommand{\arccot}{\mathrm{arccot }} % Acotanh
+
+\newcommand{\limz}[1]{\lim\limits_{#1 \rightarrow 0}} % Limes of something towards zero
+
+\newcommand{\bm}{\boldmath} % Bold font in math
+\newcommand{\dps}{\displaystyle}
+
+\newcommand{\e}{\mbox{e}} % e noncursive in math mode
+
+\newcommand{\del}{\partial} % partial diff operator
+\newcommand{\de}{\mathrm{d}} % differential d
+\newcommand{\D}{\mathrm{d}} % differential d
+\newcommand{\GRAD}{\mathrm{grad}\ } % gradient
+\newcommand{\DIV}{\mathrm{div}\ } % divergence
+\newcommand{\ROT}{\mathrm{rot}\ } % rotation
+
+\newcommand{\CONST}{\mathrm{const.\ }} % constant
+\newcommand{\var}{\mathrm{var}} % variance
+
+\newcommand{\g}{^\circ} % degrees
+\newcommand{\degr}{^\circ} % degrees
+
+\newcommand{\msol}{M_\odot} % solar mass
+\newcommand{\Lsol}{L_\odot} % solar luminosity
+\newcommand{\Lsun}{L_\odot} % solar luminosity
+\newcommand{\order}{\mathcal{O}} % order, e.g. O(h^2)
+
+
+\newcommand{\x}{\mathbf{x}} % x vector
+\newcommand{\xdot}{\dot{\mathbf{x}}} % x dot vector
+\newcommand{\xddot}{\ddot{\mathbf{x}}} % x doubledot vector
+\newcommand{\R}{\mathbf{r}} % r vector
+\newcommand{\rdot}{\dot{\mathbf{r}}} % r dot vector
+\newcommand{\rddot}{\ddot{\mathbf{r}}} % r doubledot vector
+\newcommand{\vel}{\mathbf{v}} % v vector
+\newcommand{\V}{\mathbf{v}} % v vector
+\newcommand{\vdot}{\dot{\mathbf{v}}} % v dot vector
+\newcommand{\vddot}{\ddot{\mathbf{v}}} % v doubledot vector
+
+\newcommand{\dete}{\mathrm{d}t} % dt
+\newcommand{\delte}{\del t} % partial t
+\newcommand{\dex}{\mathrm{d}x} % dx
+\newcommand{\delx}{\del x} % partial x
+\newcommand{\der}{\mathrm{d}r} % dr
+\newcommand{\delr}{\del r} % partial r
+
+
+\newcommand{\deldx}{\frac{\del}{\del x}} % shortcut partial derivative, in line
+\newcommand{\ddx}{\frac{\de}{\de x}} % shortcut total derivative, in line
+\newcommand{\DELDX}[1]{\frac{\del #1}{\del x}} % shortcut partial derivative, on fraction
+\newcommand{\DDX}[1]{\frac{\de #1}{\de x}} % shortcut total derivative, on fraction
+
+\newcommand{\deldvecx}{\frac{\del}{\del \x}} % shortcut partial derivative, in line
+\newcommand{\ddvecx}{\frac{\de}{\de \x}} % shortcut total derivative, in line
+\newcommand{\DELDVECX}[1]{\frac{\del #1}{\del \x}} % shortcut partial derivative, on fraction
+\newcommand{\DDVECX}[1]{\frac{\de #1}{\de \x}} % shortcut total derivative, on fraction
+
+\newcommand{\deldr}{\frac{\del}{\del r}} % shortcut partial derivative, in line
+\newcommand{\ddr}{\frac{\de}{\de r}} % shortcut total derivative, in line
+\newcommand{\DELDR}[1]{\frac{\del #1}{\del r}} % shortcut partial derivative, on fraction
+\newcommand{\DDR}[1]{\frac{\de #1}{\de r}} % shortcut total derivative, on fraction
+
+\newcommand{\deldt}{\frac{\del}{\del t}} % shortcut partial derivative, in line
+\newcommand{\ddt}{\frac{\de}{\de t}} % shortcut total derivative, in line
+\newcommand{\DELDT}[1]{\frac{\del #1}{\del t}} % shortcut partial derivative, on fraction
+\newcommand{\DDT}[1]{\frac{\de #1}{\de t}} % shortcut total derivative, on fraction
+
+\newcommand{\deldxalpha}{\frac{\del}{\del x^\alpha}}
+\newcommand{\DELDXALPHA}[1]{\frac{\del #1}{\del x^\alpha}}
+\newcommand{\ddxalpha}{\frac{\de}{\de x^\alpha}}
+\newcommand{\DDXALPHA}[1]{\frac{\de #1}{\de x^\alpha}}
+
+
+
+
+
diff --git a/theory/RadiativeTransfer/GEARRT/my_defines.tex b/theory/RadiativeTransfer/GEARRT/my_defines.tex
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/theory/RadiativeTransfer/GEARRT/references.bib b/theory/RadiativeTransfer/GEARRT/references.bib
new file mode 100644
index 0000000000000000000000000000000000000000..957f446b944c2c6506e065aeeb60de4e615a9d1f
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/references.bib
@@ -0,0 +1,430 @@
+
+@article{hopkinsGIZMONewClass2015,
+ archivePrefix = {arXiv},
+ eprinttype = {arxiv},
+ eprint = {1409.7395},
+ title = {{{GIZMO}}: {{A New Class}} of {{Accurate}}, {{Mesh}}-{{Free Hydrodynamic Simulation Methods}}},
+ volume = {450},
+ issn = {0035-8711, 1365-2966},
+ shorttitle = {{{GIZMO}}},
+ abstract = {We present two new Lagrangian methods for hydrodynamics, in a systematic comparison with moving-mesh, SPH, and stationary (non-moving) grid methods. The new methods are designed to simultaneously capture advantages of both smoothed-particle hydrodynamics (SPH) and grid-based/adaptive mesh refinement (AMR) schemes. They are based on a kernel discretization of the volume coupled to a high-order matrix gradient estimator and a Riemann solver acting over the volume 'overlap.' We implement and test a parallel, second-order version of the method with self-gravity \& cosmological integration, in the code GIZMO: this maintains exact mass, energy and momentum conservation; exhibits superior angular momentum conservation compared to all other methods we study; does not require 'artificial diffusion' terms; and allows the fluid elements to move with the flow so resolution is automatically adaptive. We consider a large suite of test problems, and find that on all problems the new methods appear competitive with moving-mesh schemes, with some advantages (particularly in angular momentum conservation), at the cost of enhanced noise. The new methods have many advantages vs. SPH: proper convergence, good capturing of fluid-mixing instabilities, dramatically reduced 'particle noise' \& numerical viscosity, more accurate sub-sonic flow evolution, \& sharp shock-capturing. Advantages vs. non-moving meshes include: automatic adaptivity, dramatically reduced advection errors \& numerical overmixing, velocity-independent errors, accurate coupling to gravity, good angular momentum conservation and elimination of 'grid alignment' effects. We can, for example, follow hundreds of orbits of gaseous disks, while AMR and SPH methods break down in a few orbits. However, fixed meshes minimize 'grid noise.' These differences are important for a range of astrophysical problems.},
+ number = {1},
+ journal = {Monthly Notices of the Royal Astronomical Society},
+ doi = {10.1093/mnras/stv195},
+ author = {Hopkins, Philip F.},
+ month = jun,
+ year = {2015},
+ keywords = {Astrophysics - Cosmology and Nongalactic Astrophysics,Astrophysics - Instrumentation and Methods for Astrophysics,Physics - Computational Physics,Astrophysics - Astrophysics of Galaxies,Physics - Fluid Dynamics,_tablet},
+ pages = {53-110},
+ file = {/home/mivkov/Zotero/storage/23ALA3M9/Hopkins_2015_GIZMO.pdf;/home/mivkov/Zotero/storage/SHV6QDMB/1409.html}
+}
+
+@article{lansonRenormalizedMeshfreeSchemes2008,
+ title = {Renormalized {{Meshfree Schemes I}}: {{Consistency}}, {{Stability}}, and {{Hybrid Methods}} for {{Conservation Laws}}},
+ volume = {46},
+ issn = {0036-1429},
+ shorttitle = {Renormalized {{Meshfree Schemes I}}},
+ abstract = {This paper is devoted to the study of a new kind of meshfree scheme based on a new class of meshfree derivatives: the renormalized meshfree derivatives, which improve the consistency of the original weighted particle methods. The weak renormalized meshfree scheme, built from the weak formulation of general conservation laws, turns out to be \$L\^2\$ stable under some geometrical conditions on the distribution of particles and some regularity conditions of the transport field. A time discretization is then performed by analogy with finite volume methods, and the \$L\^1\$, \$L\^\textbackslash{}infty\$, and \$BV\$ stabilities of the obtained time discretized scheme are studied. From the same analogy with finite volume methods, a hybrid particle scheme is built using the Godunov method and is numerically compared to the weak renormalized scheme.},
+ number = {4},
+ journal = {SIAM J. Numer. Anal.},
+ doi = {10.1137/S0036142903427718},
+ author = {Lanson, Nathalie and Vila, Jean-Paul},
+ month = apr,
+ year = {2008},
+ keywords = {finite volume,hybrid particle schemes,meshfree methods,nonlinear conservation law,renormalization,SPH,weak scheme,_tablet},
+ pages = {1912--1934},
+ file = {/home/mivkov/Zotero/storage/QPBH9IJ6/Lanson_Vila_2008_Renormalized Meshfree Schemes I.pdf}
+}
+
+@article{ivanovaCommonEnvelopeEvolution2013,
+ archivePrefix = {arXiv},
+ eprinttype = {arxiv},
+ eprint = {1209.4302},
+ title = {Common {{Envelope Evolution}}: {{Where}} We Stand and How We Can Move Forward},
+ volume = {21},
+ issn = {0935-4956, 1432-0754},
+ shorttitle = {Common {{Envelope Evolution}}},
+ abstract = {This work aims to present our current best physical understanding of common-envelope evolution (CEE). We highlight areas of consensus and disagreement, and stress ideas which should point the way forward for progress in this important but long-standing and largely unconquered problem. Unusually for CEE-related work, we mostly try to avoid relying on results from population synthesis or observations, in order to avoid potentially being misled by previous misunderstandings. As far as possible we debate all the relevant issues starting from physics alone, all the way from the evolution of the binary system immediately before CEE begins to the processes which might occur just after the ejection of the envelope. In particular, we include extensive discussion about the energy sources and sinks operating in CEE, and hence examine the foundations of the standard energy formalism. Special attention is also given to comparing the results of hydrodynamic simulations from different groups and to discussing the potential effect of initial conditions on the differences in the outcomes. We compare current numerical techniques for the problem of CEE and also whether more appropriate tools could and should be produced (including new formulations of computational hydrodynamics, and attempts to include 3D processes within 1D codes). Finally we explore new ways to link CEE with observations. We compare previous simulations of CEE to the recent outburst from V1309 Sco, and discuss to what extent post-common-envelope binaries and nebulae can provide information, e.g. from binary eccentricities, which is not currently being fully exploited.},
+ number = {1},
+ journal = {The Astronomy and Astrophysics Review},
+ doi = {10.1007/s00159-013-0059-2},
+ author = {Ivanova, N. and Justham, S. and Chen, X. and De Marco, O. and Fryer, C. L. and Gaburov, E. and Ge, H. and Glebbeek, E. and Han, Z. and Li, X.-D. and Lu, G. and Marsh, T. and Podsiadlowski, Ph and Potter, A. and Soker, N. and Taam, R. and Tauris, T. M. and van den Heuvel, E. P. J. and Webbink, R. F.},
+ month = nov,
+ year = {2013},
+ keywords = {Astrophysics - High Energy Astrophysical Phenomena,Astrophysics - Solar and Stellar Astrophysics,_tablet},
+ file = {/home/mivkov/Zotero/storage/V6IL3NUY/Ivanova et al_2013_Common Envelope Evolution.pdf;/home/mivkov/Zotero/storage/MMFFSMV7/1209.html}
+}
+
+@phdthesis{vandenbrouckeAdvancedModelsSimulating,
+ title = {{Advanced models for simulating dwarf galaxy formation and evolution}},
+ language = {nl},
+ author = {Vandenbroucke, Bert},
+ year = {2016},
+ file = {/home/mivkov/Zotero/storage/T3UZJAWQ/Vandenbroucke - Advanced models for simulating dwarf galaxy format.pdf}
+}
+
+
+
+@article{ramses-rt13,
+ title = {{{RAMSES}}-{{RT}}: Radiation Hydrodynamics in the Cosmological Context},
+ shorttitle = {{{RAMSES}}-{{RT}}},
+ author = {Rosdahl, J. and Blaizot, J. and Aubert, D. and Stranex, T. and Teyssier, R.},
+ year = {2013},
+ month = dec,
+ volume = {436},
+ pages = {2188--2231},
+ issn = {0035-8711},
+ doi = {10.1093/mnras/stt1722},
+ abstract = {We present a new implementation of radiation hydrodynamics (RHD) in the adaptive mesh refinement (AMR) code RAMSES. The multigroup radiative transfer (RT) is performed on the AMR grid with a first-order Godunov method using the M1 closure for the Eddington tensor, and is coupled to the hydrodynamics via non-equilibrium thermochemistry of hydrogen and helium. This moment-based approach has the great advantage that the computational cost is independent of the number of radiative sources - it can even deal with continuous regions of emission such as bound-free emission from gas. As it is built directly into RAMSES, the RT takes natural advantage of the refinement and parallelization strategies already in place. Since we use an explicit advection solver for the radiative transport, the time-step is restricted by the speed of light - a severe limitation that can be alleviated using the so-called reduced speed of light approximation. We propose a rigorous framework to assess the validity of this approximation in various conditions encountered in cosmology and galaxy formation. We finally perform with our newly developed code a complete suite of RHD tests, comparing our results to other RHD codes. The tests demonstrate that our code performs very well and is ideally suited for exploring the effect of radiation on current scenarios of structure and galaxy formation.},
+ file = {/home/mivkov/Zotero/storage/M5QSAVLR/Rosdahl et al_2013_RAMSES-RT.pdf},
+ journal = {Monthly Notices of the Royal Astronomical Society},
+ keywords = {methods: numerical,radiative transfer}
+}
+
+@article{ramses-rt15,
+ title = {A Scheme for Radiation Pressure and Photon Diffusion with the {{M1}} Closure in {{RAMSES}}-{{RT}}},
+ author = {Rosdahl, J. and Teyssier, R.},
+ year = {2015},
+ month = jun,
+ volume = {449},
+ pages = {4380--4403},
+ issn = {0035-8711, 1365-2966},
+ doi = {10.1093/mnras/stv567},
+ abstract = {We describe and test an updated version of radiation-hydrodynamics (RHD) in the RAMSES code, that includes three new features: i) radiation pressure on gas, ii) accurate treatment of radiation diffusion in an unresolved optically thick medium, and iii) relativistic corrections that account for Doppler effects and work done by the radiation to first order in v/c. We validate the implementation in a series of tests, which include a morphological assessment of the M1 closure for the Eddington tensor in an astronomically relevant setting, dust absorption in a optically semi-thick medium, direct pressure on gas from ionising radiation, convergence of our radiation diffusion scheme towards resolved optical depths, correct diffusion of a radiation flash and a constant luminosity radiation, and finally, an experiment from Davis et al. of the competition between gravity and radiation pressure in a dusty atmosphere, and the formation of radiative Rayleigh-Taylor instabilities. With the new features, RAMSES-RT can be used for state-of-the-art simulations of radiation feedback from first principles, on galactic and cosmological scales, including not only direct radiation pressure from ionising photons, but also indirect pressure via dust from multi-scattered IR photons reprocessed from higher-energy radiation, both in the optically thin and thick limits.},
+ archivePrefix = {arXiv},
+ eprint = {1411.6440},
+ eprinttype = {arxiv},
+ file = {/home/mivkov/Zotero/storage/8WS8SH75/Rosdahl_Teyssier_2015_A scheme for radiation pressure and photon diffusion with the M1 closure in.pdf;/home/mivkov/Zotero/storage/U86AN9F6/1411.html},
+ journal = {Monthly Notices of the Royal Astronomical Society},
+ keywords = {_tablet_modified,Astrophysics - High Energy Astrophysical Phenomena,Astrophysics - Instrumentation and Methods for Astrophysics},
+ number = {4}
+}
+
+
+@article{ilievCosmologicalRadiativeTransfer2006,
+ title = {Cosmological {{Radiative Transfer Codes Comparison Project I}}: {{The Static Density Field Tests}}},
+ shorttitle = {Cosmological {{Radiative Transfer Codes Comparison Project I}}},
+ author = {Iliev, Ilian T. and Ciardi, Benedetta and Alvarez, Marcelo A. and Maselli, Antonella and Ferrara, Andrea and Gnedin, Nickolay Y. and Mellema, Garrelt and Nakamoto, Taishi and Norman, Michael L. and Razoumov, Alexei O. and Rijkhorst, Erik-Jan and Ritzerveld, Jelle and Shapiro, Paul R. and Susa, Hajime and Umemura, Masayuki and Whalen, Daniel J.},
+ year = {2006},
+ month = mar,
+ journal = {arXiv:astro-ph/0603199},
+ eprint = {astro-ph/0603199},
+ eprinttype = {arxiv},
+ doi = {10.1111/j.1365-2966.2006.10775.x},
+ abstract = {Radiative transfer simulations are now at the forefront of numerical astrophysics. They are becoming crucial for an increasing number of astrophysical and cosmological problems; at the same time their computational cost has come to the reach of currently available computational power. Further progress is retarded by the considerable number of different algorithms (including various flavours of ray-tracing and moment schemes) developed, which makes the selection of the most suitable technique for a given problem a non-trivial task. Assessing the validity ranges, accuracy and performances of these schemes is the main aim of this paper, for which we have compared 11 independent RT codes on 5 test problems: (0) basic physics, (1) isothermal H II region expansion and (2) H II region expansion with evolving temperature, (3) I-front trapping and shadowing by a dense clump, (4) multiple sources in a cosmological density field. The outputs of these tests have been compared and differences analyzed. The agreement between the various codes is satisfactory although not perfect. The main source of discrepancy appears to reside in the multi-frequency treatment approach, resulting in different thicknesses of the ionized-neutral transition regions and the temperature structure. The present results and tests represent the most complete benchmark available for the development of new codes and improvement of existing ones. To this aim all test inputs and outputs are made publicly available in digital form.},
+ archiveprefix = {arXiv},
+ langid = {english},
+ keywords = {Astrophysics},
+ file = {/home/mivkov/Zotero/storage/6L4AAYAN/Iliev et al. - 2006 - Cosmological Radiative Transfer Codes Comparison P.pdf}
+}
+
+
+
+@article{gonzalezHERACLESThreedimensionalRadiation2007,
+ title = {{{HERACLES}}: A Three-Dimensional Radiation Hydrodynamics Code},
+ shorttitle = {{{HERACLES}}},
+ author = {Gonz{\'a}lez, M. and Audit, E. and Huynh, P.},
+ year = {2007},
+ month = mar,
+ journal = {A\&A},
+ volume = {464},
+ number = {2},
+ pages = {429--435},
+ issn = {0004-6361, 1432-0746},
+ doi = {10.1051/0004-6361:20065486},
+ abstract = {Methods. The radiation transfer is modelled using a two-moment model and a closure
+relation that allows large angular anisotropies in the radiation field to be preserved and
+reproduced. The radiative equations thus obtained are solved by a second-order Godunov-type method
+and integrated implicitly by using iterative solvers. HERACLES has been parallelized with the MPI
+library and implemented in Cartesian, cylindrical, and spherical coordinates. To characterize the
+accuracy of HERACLES and to compare it with other codes, we performed a series of tests including
+purely radiative tests and radiation-hydrodynamics ones. Results. The results show that the physical
+model used in HERACLES for the transfer is fairly accurate in both the diffusion and transport
+limit, but also for semi-transparent regions. Conclusions. This makes HERACLES very well-suited to
+studying many astrophysical problems such as radiative shocks, molecular jets of young stars,
+fragmentation and formation of dense cores in the interstellar medium, and protoplanetary discs.},
+ langid = {english},
+ file = {/home/mivkov/Zotero/storage/NKTGBME2/González et al. - 2007 - HERACLES a three-dimensional
+radiation hydrodynam.pdf}
+}
+
+@article{ilievCosmologicalRadiativeTransfer2006,
+ title = {Cosmological {{Radiative Transfer Codes Comparison Project I}}: {{The Static Density
+Field Tests}}},
+ shorttitle = {Cosmological {{Radiative Transfer Codes Comparison Project I}}},
+ author = {Iliev, Ilian T. and Ciardi, Benedetta and Alvarez, Marcelo A. and Maselli, Antonella and
+Ferrara, Andrea and Gnedin, Nickolay Y. and Mellema, Garrelt and Nakamoto, Taishi and Norman,
+Michael L. and Razoumov, Alexei O. and Rijkhorst, Erik-Jan and Ritzerveld, Jelle and Shapiro, Paul
+R. and Susa, Hajime and Umemura, Masayuki and Whalen, Daniel J.},
+ year = {2006},
+ month = mar,
+ journal = {arXiv:astro-ph/0603199},
+ eprint = {astro-ph/0603199},
+ eprinttype = {arxiv},
+ doi = {10.1111/j.1365-2966.2006.10775.x},
+ abstract = {Radiative transfer simulations are now at the forefront of numerical astrophysics.
+They are becoming crucial for an increasing number of astrophysical and cosmological problems; at
+the same time their computational cost has come to the reach of currently available computational
+power. Further progress is retarded by the considerable number of different algorithms (including
+various flavours of ray-tracing and moment schemes) developed, which makes the selection of the most
+suitable technique for a given problem a non-trivial task. Assessing the validity ranges, accuracy
+and performances of these schemes is the main aim of this paper, for which we have compared 11
+independent RT codes on 5 test problems: (0) basic physics, (1) isothermal H II region expansion and
+(2) H II region expansion with evolving temperature, (3) I-front trapping and shadowing by a dense
+clump, (4) multiple sources in a cosmological density field. The outputs of these tests have been
+compared and differences analyzed. The agreement between the various codes is satisfactory although
+not perfect. The main source of discrepancy appears to reside in the multi-frequency treatment
+approach, resulting in different thicknesses of the ionized-neutral transition regions and the
+temperature structure. The present results and tests represent the most complete benchmark available
+for the development of new codes and improvement of existing ones. To this aim all test inputs and
+outputs are made publicly available in digital form.},
+ archiveprefix = {arXiv},
+ langid = {english},
+ keywords = {Astrophysics},
+ file = {/home/mivkov/Zotero/storage/6L4AAYAN/Iliev et al. - 2006 - Cosmological Radiative Transfer
+Codes Comparison P.pdf}
+}
+
+@article{ilievCosmologicalRadiativeTransfer2009,
+ title = {Cosmological {{Radiative Transfer Comparison Project II}}: {{The Radiation-Hydrodynamic
+Tests}}},
+ shorttitle = {Cosmological {{Radiative Transfer Comparison Project II}}},
+ author = {Iliev, Ilian T. and Whalen, Daniel and Mellema, Garrelt and Ahn, Kyungjin and Baek,
+Sunghye and Gnedin, Nickolay Y. and Kravtsov, Andrey V. and Norman, Michael and Raicevic, Milan and
+Reynolds, Daniel R. and Sato, Daisuke and Shapiro, Paul R. and Semelin, Benoit and Smidt, Joseph and
+Susa, Hajime and Theuns, Tom and Umemura, Masayuki},
+ year = {2009},
+ month = dec,
+ journal = {Monthly Notices of the Royal Astronomical Society},
+ volume = {400},
+ number = {3},
+ eprint = {0905.2920},
+ eprinttype = {arxiv},
+ pages = {1283--1316},
+ issn = {00358711, 13652966},
+ doi = {10.1111/j.1365-2966.2009.15558.x},
+ abstract = {The development of radiation hydrodynamical methods that are able to follow gas
+dynamics and radiative transfer self-consistently is key to the solution of many problems in
+numerical astrophysics. Such fluid flows are highly complex, rarely allowing even for approximate
+analytical solutions against which numerical codes can be tested. An alternative validation
+procedure is to compare different methods against each other on common problems, in order to assess
+the robustness of the results and establish a range of validity for the methods. Previously, we
+presented such a comparison for a set of pure radiative transfer tests (i.e. for fixed, non-evolving
+density fields). This is the second paper of the Cosmological Radiative Transfer (RT) Comparison
+Project, in which we compare 9 independent RT codes directly coupled to gasdynamics on 3 relatively
+simple astrophysical hydrodynamics problems: (5) the expansion of an H II region in a uniform
+medium; (6) an ionization front (I-front) in a 1/r2 density profile with a flat core, and (7), the
+photoevaporation of a uniform dense clump. Results show a broad agreement between the different
+methods and no big failures, indicating that the participating codes have reached a certain level of
+maturity and reliability. However, many details still do differ, and virtually every code has showed
+some shortcomings and has disagreed, in one respect or another, with the majority of the results.
+This underscores the fact that no method is universal and all require careful testing of the
+particular features which are most relevant to the specific problem at hand.},
+ archiveprefix = {arXiv},
+ langid = {english},
+ keywords = {Astrophysics - Cosmology and Nongalactic Astrophysics},
+ file = {/home/mivkov/Zotero/storage/2YNJDDST/Iliev et al. - 2009 - Cosmological Radiative Transfer
+Comparison Project.pdf}
+}
+
+@article{katzCosmologicalSimulationsTreeSPH1996,
+ title = {Cosmological {{Simulations}} with {{TreeSPH}}},
+ author = {Katz, Neal and Weinberg, David H. and Hernquist, Lars},
+ year = {1996},
+ month = jul,
+ journal = {The Astrophysical Journal Supplement Series},
+ volume = {105},
+ pages = {19},
+ issn = {0067-0049},
+ doi = {10.1086/192305},
+ abstract = {We describe numerical methods for incorporating gasdynamics into cosmological
+simulations and present illustrative applications to the cold dark matter (CDM) scenario. Our
+evolution code, a version of TreeSPH (Hernquist \& Katz 1989) generalized to handle comoving
+coordinates and periodic boundary conditions, combines smoothed-particle hydrodynamics (SPH) with
+the hierarchical tree method for computing gravitational forces. The Lagrangian hydrodynamics
+approach and individual time steps for gas particles give the algorithm a large dynamic range, which
+is essential for studies of galaxy formation in a cosmological context. The code incorporates
+radiative cooling for an optically thin, primordial composition gas in ionization equilibrium with a
+user-specified ultraviolet background. We adopt a phenomenological prescription for star formation
+that gradually turns cold, dense, Jeans-unstable gas into collisionless stars, returning supernova
+feedback energy to the surrounding medium. In CDM simulations, some of the baryons that fall into
+dark matter potential wells dissipate their acquired thermal energy and condense into clumps with
+roughly galactic masses. The resulting galaxy population is insensitive to assumptions about star
+formation; we obtain similar baryonic mass functions and galaxy correlation functions from
+simulations with star formation and from simulations without star formation in which we identify
+galaxies directly from the cold, dense gas.},
+ keywords = {_tablet,COSMOLOGY: DARK MATTER,COSMOLOGY: LARGE-SCALE STRUCTURE OF UNIVERSE,COSMOLOGY:
+THEORY,GALAXIES: FORMATION,HYDRODYNAMICS,METHODS: NUMERICAL},
+ file = {/home/mivkov/Zotero/storage/2YV7DC26/Katz et al_1996_Cosmological Simulations with
+TreeSPH.pdf}
+}
+
+@article{levermoreRelatingEddingtonFactors1984a,
+ title = {Relating {{Eddington}} Factors to Flux Limiters},
+ author = {Levermore, C. D.},
+ year = {1984},
+ month = feb,
+ journal = {Journal of Quantitative Spectroscopy and Radiative Transfer},
+ volume = {31},
+ number = {2},
+ pages = {149--160},
+ issn = {0022-4073},
+ doi = {10.1016/0022-4073(84)90112-2},
+ abstract = {Variable Eddington factors and flux-limiters have been introduced in the P-1 and
+diffusion equations, respectively, to handle situations when the underlying intensity is too
+anisotropic for the unmodified theories to remain valid. We present a derivation of a relation
+between the two for which a new approach to the diffusion approximation is used. Algebraic
+expressions for Eddington factors satisfying the moment conditions are not satisfactory for closing
+the P-1 equations but, by using the derived relation, yield acceptable flux-limited diffusion
+theories.},
+ langid = {english},
+ file = {/home/mivkov/Zotero/storage/WSX9USJ7/Levermore_1984_Relating Eddington factors to flux
+limiters.pdf;/home/mivkov/Zotero/storage/XAJD9TFY/0022407384901122.html}
+}
+
+@article{smithGrackleChemistryCooling2017a,
+ title = {Grackle: A {{Chemistry}} and {{Cooling Library}} for {{Astrophysics}}},
+ shorttitle = {Grackle},
+ author = {Smith, Britton D. and Bryan, Greg L. and Glover, Simon C. O. and Goldbaum, Nathan J. and
+Turk, Matthew J. and Regan, John and Wise, John H. and Schive, Hsi-Yu and Abel, Tom and Emerick,
+Andrew and O'Shea, Brian W. and Anninos, Peter and Hummels, Cameron B. and Khochfar, Sadegh},
+ year = {2017},
+ month = apr,
+ journal = {Mon. Not. R. Astron. Soc.},
+ volume = {466},
+ number = {2},
+ eprint = {1610.09591},
+ eprinttype = {arxiv},
+ pages = {2217--2234},
+ issn = {0035-8711, 1365-2966},
+ doi = {10.1093/mnras/stw3291},
+ abstract = {We present the Grackle chemistry and cooling library for astrophysical simulations and
+models. Grackle provides a treatment of non-equilibrium primordial chemistry and cooling for H, D,
+and He species, including H2 formation on dust grains; tabulated primordial and metal cooling;
+multiple UV background models; and support for radiation transfer and arbitrary heat sources. The
+library has an easily implementable interface for simulation codes written in C, C++, and Fortran as
+well as a Python interface with added convenience functions for semi-analytical models. As an
+open-source project, Grackle provides a community resource for accessing and disseminating
+astrochemical data and numerical methods. We present the full details of the core functionality, the
+simulation and Python interfaces, testing infrastructure, performance, and range of applicability.
+Grackle is a fully open-source project and new contributions are welcome.},
+ archiveprefix = {arXiv},
+ keywords = {Astrophysics - Astrophysics of Galaxies,Astrophysics - Cosmology and Nongalactic
+Astrophysics,Astrophysics - Instrumentation and Methods for Astrophysics},
+ file = {/home/mivkov/Zotero/storage/8YEYQAV2/Smith et
+al_2017_Grackle.pdf;/home/mivkov/Zotero/storage/QRXIKLGA/1610.html}
+}
+
+
+@article{ivkovicThesis,
+ author = {{Ivkovic}, Mladen},
+ title = "{GEAR-RT: Towards Exa-Scale Moment Based Radiative Transfer For Cosmological
+Simulations Using Task-Based Parallelism And Dynamic Sub-Cycling with SWIFT}",
+ journal = {arXiv e-prints},
+ keywords = {Astrophysics - Instrumentation and Methods for Astrophysics},
+ year = 2023,
+ month = feb,
+ eid = {arXiv:2302.12727},
+ pages = {arXiv:2302.12727},
+ doi = {10.48550/arXiv.2302.12727},
+archivePrefix = {arXiv},
+ eprint = {2302.12727},
+ primaryClass = {astro-ph.IM},
+ adsurl = {https://ui.adsabs.harvard.edu/abs/2023arXiv230212727I},
+ adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+
+@article{vernerAtomicDataAstrophysics1996,
+ title = {Atomic {{Data}} for {{Astrophysics}}. {{II}}. {{New Analytic Fits}} for {{Photoionization
+Cross Sections}} of {{Atoms}} and {{Ions}}},
+ author = {Verner, D. A. and Ferland, G. J. and Korista, K. T. and Yakovlev, D. G.},
+ year = {1996},
+ month = jul,
+ journal = {The Astrophysical Journal},
+ volume = {465},
+ eprint = {astro-ph/9601009},
+ pages = {487},
+ issn = {0004-637X, 1538-4357},
+ doi = {10.1086/177435},
+ urldate = {2022-11-11},
+ abstract = {We present a complete set of analytic fits to the non-relativistic photoionization
+cross sections for the ground states of atoms and ions of elements from H through Si, and S, Ar, Ca,
+and Fe. Near the ionization thresholds, the fits are based on the Opacity Project theoretical cross
+sections interpolated and smoothed over resonances. At higher energies, the fits reproduce
+calculated Hartree-Dirac-Slater photoionization cross sections.},
+ archiveprefix = {arxiv},
+ langid = {english},
+ keywords = {Astrophysics,Physics - Atomic Physics},
+ file = {/home/mivkov/Zotero/storage/GAQCFADD/Verner et al. - 1996 - Atomic Data for Astrophysics.
+II. New Analytic Fit.pdf}
+}
+
+
+
+@article{smithGRACKLEChemistryCooling2017,
+ ids = {smithGRACKLEChemistryCooling2017a,smithGrackleChemistryCooling2017a},
+ title = {{{GRACKLE}}: A Chemistry and Cooling Library for Astrophysics},
+ shorttitle = {{{GRACKLE}}},
+ author = {Smith, Britton D. and Bryan, Greg L. and Glover, Simon C. O. and Goldbaum, Nathan J. and
+Turk, Matthew J. and Regan, John and Wise, John H. and Schive, Hsi-Yu and Abel, Tom and Emerick,
+Andrew and O'Shea, Brian W. and Anninos, Peter and Hummels, Cameron B. and Khochfar, Sadegh},
+ year = {2017},
+ month = apr,
+ journal = {Monthly Notices of the Royal Astronomical Society},
+ volume = {466},
+ eprint = {1610.09591},
+ pages = {2217--2234},
+ issn = {0035-8711},
+ doi = {10.1093/mnras/stw3291},
+ urldate = {2021-10-27},
+ abstract = {We present the GRACKLE chemistry and cooling library for astrophysical simulations and
+models. GRACKLE provides a treatment of non-equilibrium primordial chemistry and cooling for H, D
+and He species, including H2 formation on dust grains; tabulated primordial and metal cooling;
+multiple ultraviolet background models; and support for radiation transfer and arbitrary heat
+sources. The library has an easily implementable interface for simulation codes written in C, C++
+and FORTRAN as well as a PYTHON interface with added convenience functions for semi-analytical
+models. As an open-source project, GRACKLE provides a community resource for accessing and
+disseminating astrochemical data and numerical methods. We present the full details of the core
+functionality, the simulation and PYTHON interfaces, testing infrastructure, performance and range
+of applicability. GRACKLE is a fully open-source project and new contributions are welcome.},
+ archiveprefix = {arxiv},
+ keywords = {astrochemistry,Astrophysics - Astrophysics of Galaxies,Astrophysics - Cosmology and
+Nongalactic Astrophysics,Astrophysics - Instrumentation and Methods for Astrophysics,galaxies:
+formation,methods: numerical},
+ annotation = {ADS Bibcode: 2017MNRAS.466.2217S},
+ file = {/home/mivkov/Zotero/storage/7JY7484L/Smith et
+al_2017_GRACKLE.pdf;/home/mivkov/Zotero/storage/5MWMNWRG/abstract.html;/home/mivkov/Zotero/storage/
+QRXIKLGA/1610.html}
+}
+
+
+@article{gnedinMultidimensionalCosmologicalRadiative2001,
+ title = {Multi-Dimensional Cosmological Radiative Transfer with a {{Variable Eddington
+Tensor}} Formalism},
+ author = {Gnedin, Nickolay Y. and Abel, Tom},
+ year = {2001},
+ month = oct,
+ journal = {New Astronomy},
+ volume = {6},
+ pages = {437--455},
+ issn = {1384-1076},
+ doi = {10.1016/S1384-1076(01)00068-9},
+ urldate = {2022-11-14},
+ abstract = {We present a new approach to numerically model continuum radiative transfer based on
+the Optically Thin Variable Eddington Tensor (OTVET) approximation. Our method insures the exact
+conservation of the photon number and flux (in the explicit formulation) and automatically switches
+from the optically thick to the optically thin regime. It scales as N log N with the number of
+hydrodynamic resolution elements and is independent of the number of sources of ionizing radiation
+(i.e. works equally fast for an arbitrary source function). We also describe an implementation of
+the algorithm in a Soften Lagrangian Hydrodynamic code (SLH) and a multi-frequency approach
+appropriate for hydrogen and helium continuum opacities. We present extensive tests of our method
+for single and multiple sources in homogeneous and inhomogeneous density distributions, as well as a
+realistic simulation of cosmological reionization.},
+ keywords = {Astrophysics},
+ annotation = {ADS Bibcode: 2001NewA....6..437G},
+ file = {/home/mivkov/Zotero/storage/FN64KVYB/Gnedin_Abel_2001_Multi-dimensional cosmological
+radiative transfer with a Variable Eddington.pdf}
+}
+
+
diff --git a/theory/RadiativeTransfer/GEARRT/run.sh b/theory/RadiativeTransfer/GEARRT/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..11be530f1e14df5686245a9fa5a07e8bf94d42eb
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/run.sh
@@ -0,0 +1,7 @@
+#!/bin/bash
+echo "Generating PDF..."
+
+pdflatex -jobname=GEARRT GEARRT.tex
+bibtex GEARRT
+pdflatex -jobname=GEARRT GEARRT.tex
+pdflatex -jobname=GEARRT GEARRT.tex
diff --git a/theory/RadiativeTransfer/GEARRT/tables/fitting_parameters.tex b/theory/RadiativeTransfer/GEARRT/tables/fitting_parameters.tex
new file mode 100644
index 0000000000000000000000000000000000000000..48aee571b4ba0c8bd3716444d86d35d281e2ee3e
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/tables/fitting_parameters.tex
@@ -0,0 +1,16 @@
+
+\begin{table}
+\centering
+\begin{tabular}{l | c c c c c c c }
+ & $E_0$ [eV] & $\sigma_0$ [cm$^2$] & $y_a$ & $P$ & $y_w$ & $y_0$ & $y_1$ \\
+\hline\\[-0.5em]
+H$^0$ & 0.4298 & 5.475$\times 10^{-14}$ & 32.88 & 2.963 & 0 & 0 & 0 \\
+He$^0$ & 13.61 & 9.492$\times 10^{-16}$ & 1.469 & 3.188 & 2.039 & 0.4434 & 2.136 \\
+He$^+$ & 1.720 & 1.369$\times 10^{-14}$ & 32.88 & 2.963 & 0 & 0 & 0 \\
+ \end{tabular}
+\caption{Fitting parameters used for the ionizing cross section parametrizations given in
+eq.~\ref{eq:sigma-parametrizaiton}. The parametrization and these parameters are taken from
+\cite{vernerAtomicDataAstrophysics1996}. Note that the table given in the appendix of
+\cite{ramses-rt13} contains a typo in $E_0$ for He$^0$, and is off by a factor of 100.}
+\label{tab:cross-sections}
+\end{table}
diff --git a/theory/RadiativeTransfer/GEARRT/tables/rt_variables.tex b/theory/RadiativeTransfer/GEARRT/tables/rt_variables.tex
new file mode 100644
index 0000000000000000000000000000000000000000..7284945af25f8ce4d406dd7126fd0e06652044ed
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/tables/rt_variables.tex
@@ -0,0 +1,125 @@
+%\begin{table}
+\begin{center}
+\begin{small}
+\begin{tabular}{p{0.32\textwidth} p{0.42\textwidth} p{0.26\textwidth}}
+
+
+variable & name & units (cgs) \\[.5em]
+%---------------------------------------------------------------------------------------------------
+\hline \\
+
+$I_\nu (\x, \mathbf{n}, t)$ &
+ specific intensity &
+ $[I_\nu] = \frac{\text{erg}}{\text{cm}^2 \text{ rad Hz s}}$
+\\[.5em]
+$u_\nu (\x, \mathbf{n}, t) = \frac{I_\nu}{c}$ &
+ radiation energy density &
+ $[u_\nu] = \frac{\text{erg}}{\text{cm}^3 \text{ rad Hz}}$
+\\[.5em]
+$E_\nu (\x, t) = \int_{4 \pi} \frac{I_\nu}{c} \de \Omega$ &
+ total energy density &
+ $[E_\nu] = \frac{\text{erg}}{\text{cm}^3 \text{ Hz}}$
+\\[.5em]
+$N_\nu (\x, t) = E_\nu / h$ &
+ photon number density &
+ $[N_\nu] = \frac{1}{\text{cm}^3 \text{ Hz}}$
+\\[.5em]
+$E_{rad} (\x, t) = \int_{0}^{\infty} E_\nu \de \nu$ &
+ total integrated energy density &
+ $[E_{rad}] = \frac{\text{erg}}{\text{cm}^3}$
+\\[.5em]
+$N_{rad} (\x, t) = E_{rad} / h$ &
+ integrated photon number density &
+ $[N_{rad}] = \frac{1}{\text{cm}^3}$
+\\[.5em]
+% $J_\nu (\x, t) = \int_{4 \pi} I_\nu \frac{\de \Omega}{4 \pi}
+% = \frac{c}{4 \pi} E_\nu$ &
+% mean radiation specific intensity &
+% $[J_\nu] = \frac{\text{erg}}{\text{cm}^2 \text{ Hz s}}$
+% \\[.5em]
+$\mathbf{F}_\nu(\x, t) = \int_{4 \pi} I_\nu \mathbf{n} \de \Omega$ &
+ radiation flux &
+ $[\mathbf{F}_\nu] = \frac{\text{erg}}{\text{cm}^2 \text{ s Hz}}$
+\\[.5em]
+$\mathbf{P}_\nu (\x, t) = \frac{\F_\nu}{c^2}$ &
+ radiation momentum density &
+ $[\mathbf{P}_\nu] = \frac{\text{erg}}{ \text{cm}^4 \text{ s}^{-1} \text{ Hz}}$
+\\[.5em]
+$\mathds{P}_\nu (\x, t) = \int_{4 \pi} \frac{I_\nu}{c} \mathbf{n} \otimes \mathbf{n} \de \Omega$ &
+ radiation pressure tensor &
+ $[\mathds{P}_\nu ] = \frac{\text{erg}}{\text{cm}^3 \text{ Hz}}$
+
+
+\\
+%---------------------------------------------------------------------------------------------------
+\hline\\
+
+
+$\eta_\nu(\x, t)$ &
+ source function (of $I_\nu$)&
+ $[ \eta_\nu ] = \frac{\text{erg}}{\text{cm}^3 \text{ rad Hz s}}$
+\\[.5em]
+$n(\x, t)$ &
+ number density &
+ $[ n ] = \text{cm}^{-3}$
+\\[.5em]
+$\mathcal{J}_\nu(T)$ &
+ (specific intensity of a) spectrum &
+ $[\mathcal{J}_\nu (T)] = \frac{\text{erg}}{\text{cm}^2 \text{ rad Hz s}}$
+\\[.5em]
+$\Gamma_{\nu, j}$ &
+ photoionization rate of species $j$ &
+ $[\Gamma] = \text{s}^{-1}$
+\\[.5em]
+$\mathcal{H}_{\nu,j}$ &
+ (photo)heating rate of species $j$ &
+ $[\mathcal{H}_{\nu,j}] = \frac{\text{erg}}{s\ cm}^{3}$
+\\[.5em]
+$\alpha_\nu = \frac{1}{\lambda_\nu} = \rho \kappa_\nu = n \sigma_\nu$ &
+ absorption coefficient &
+ $[\alpha_\nu] = \text{cm}^{-1}$
+\\[.5em]
+% $\alpha_E = \frac{\int_0^\infty \alpha_\nu E_\nu \de \nu}{\int_0^\infty E_\nu \de \nu}$ &
+% energy mean of abs. coeff. &
+% $[\alpha_E] = \text{cm}^{-1}$
+% \\[.5em]
+% $\alpha_F = \frac{\int_0^\infty \alpha_\nu F_\nu \de \nu}{\int_0^\infty F_\nu \de \nu}$ &
+% flux mean of abs. coeff. &
+% $[\alpha_F] = \text{cm}^{-1}$
+% \\[.5em]
+% $j_\nu $ &
+% emission coefficient &
+% $[j_\nu] = \frac{\text{erg}}{\text{cm}^3 \text{rad Hz s}}$
+% \\[.5em]
+$\lambda_\nu $ &
+ mean free path &
+ $[\lambda_\nu] = \text{cm}$
+\\[.5em]
+$\kappa_\nu $ &
+ opacity &
+ $[\kappa_\nu] = \frac{\text{cm}^2}{\text{g}}$
+\\[.5em]
+$\tau_\nu(s) = \int_0^s \alpha_\nu(x) \de x$ &
+ optical depth &
+ $[\tau_\nu] = 1$
+\\[.5em]
+$\sigma_{j\nu}$ &
+ interacton cross section (of species $j$)&
+ $[\sigma_{j\nu}] = \text{cm}^2$
+\\[.5em]
+$\sigma_{ij}^N = \int_{\nu_{i}} N_\nu \ \sigma_{j\nu} \de \nu / \int_{\nu_i} N_\nu \de \nu$ &
+ number weighted average cross section &
+ $[\sigma_{ij}^N] = \text{cm}^2$
+\\[.5em]
+$\sigma_{ij}^E = \int_{\nu i} E_\nu \ \sigma_{j\nu} \de \nu / \int_{\nu i} E_\nu \de \nu$ &
+ energy weighted average cross section &
+ $[\sigma_{ij}^E] = \text{cm}^2$
+%
+
+\end{tabular}
+\end{small}
+\end{center}
+%\caption{Common quantities appearing in the context of radiative transfer, and their units.}
+%\label{tab:rt-variables}
+%\end{table}
+
diff --git a/theory/RadiativeTransfer/GEARRT/tables/thermochemistry_rates.tex b/theory/RadiativeTransfer/GEARRT/tables/thermochemistry_rates.tex
new file mode 100644
index 0000000000000000000000000000000000000000..0b14e3a05e5e7f503ddc08fd72e78321fb4ae08f
--- /dev/null
+++ b/theory/RadiativeTransfer/GEARRT/tables/thermochemistry_rates.tex
@@ -0,0 +1,42 @@
+\begin{table}
+\begin{center}
+\begin{tabular}{lll}
+\hline\\
+$A_{H^+} $ &
+ $8.4 \times 10^{-11} T^{-1/2} T_3^{-0.2} (1 + T_6^{0.7})^{-1}$ &
+ RR for H$^{+}$
+\\[.5em]
+$A_{d} $ &
+ $1.5 \times 10^{-10} T^{-0.6353}$ &
+ dielectronic RR for He$^{+}$
+\\[.5em]
+$A_{He^+} $ &
+ $1.9 \times 10^{-3} T^{-1.5} \mathrm{e}^{-470000/T} (1 + 0.3 \mathrm{e}^{-94000/T})$ &
+ RR for He$^{+}$
+\\[.5em]
+$A_{He^{++}} $ &
+ $3.36 \times 10^{-10} T^{-1/2} T_3^{-0.2} (1 + T_6^{0.7})^{-1}$ &
+ RR for He$^{++}$
+\\[.5em]
+\hline\\
+$\Gamma_{H^{0}} $ &
+ $5.85 \times 10^{-11} T^{1/2} \mathrm{e}^{-157809.1/T} ( 1 + T_5^{1/2})^{-1}$&
+ CIR for H$^{0}$
+\\[.5em]
+$\Gamma_{He^{0}} $ &
+ $2.38 \times 10^{-11} T^{1/2} \mathrm{e}^{-285335.4.1/T} ( 1 + T_5^{1/2})^{-1}$&
+ CIR for He$^{0}$
+\\[.5em]
+$\Gamma_{He^{+}} $ &
+ $5.68 \times 10^{-12} T^{1/2} \mathrm{e}^{-631515/T} ( 1 + T_5^{1/2})^{-1}$&
+ CIR for He$^{+}$
+\\[.5em]
+\hline
+\end{tabular}
+\end{center}
+\caption{Temperature ($T$) dependent recombination rates (RR) and collisional ionization rates
+(CIR) for Hydrogen and Helium species, adapted from \citet{katzCosmologicalSimulationsTreeSPH1996}.
+All rates are in units of cm$^3$ s$^{-1}$. $T_n$ is shorthand for $T / 10^n$K.}
+\label{tab:coll-ion-rates-katz}
+\end{table}
+