diff --git a/examples/Traphic/getGlass.sh b/examples/Traphic/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..d5c5f590ac37c9c9431d626a2ea61b0c12c1513c
--- /dev/null
+++ b/examples/Traphic/getGlass.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_64.hdf5
diff --git a/examples/Traphic/makeIC.py b/examples/Traphic/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..c967cd4a428549110539b78bce19134979d7ce5e
--- /dev/null
+++ b/examples/Traphic/makeIC.py
@@ -0,0 +1,78 @@
+################################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+################################################################################
+
+import h5py
+from numpy import *
+
+# Generates a swift IC file for the Traphic example in a periodic box
+
+fileName = "traphic.hdf5"
+
+#---------------------------------------------------
+glass = h5py.File("glassCube_64.hdf5", "r")
+
+pos = glass["/PartType0/Coordinates"][:,:]
+h = glass["/PartType0/SmoothingLength"][:]
+
+# Merge things
+numPart = size(h)
+
+# Generate extra arrays
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+m = zeros(numPart)
+u = zeros(numPart)
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [1., 1., 1.]
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+grp.attrs["Dimension"] = 3
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+#Particle group
+grp = file.create_group("/PartType0")
+grp.create_dataset('Coordinates', data=pos, dtype='d')
+grp.create_dataset('Velocities', data=v, dtype='f')
+grp.create_dataset('Masses', data=m, dtype='f')
+grp.create_dataset('SmoothingLength', data=h, dtype='f')
+grp.create_dataset('InternalEnergy', data=u, dtype='f')
+grp.create_dataset('ParticleIDs', data=ids, dtype='L')
+
+file.close()
diff --git a/examples/Traphic/run.sh b/examples/Traphic/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..5d82a5862824bdf40c2f748bb39806b140ae6cc7
--- /dev/null
+++ b/examples/Traphic/run.sh
@@ -0,0 +1,16 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e glassCube_64.hdf5 ]
+then
+ echo "Fetching initial glass file for the Traphic example..."
+ ./getGlass.sh
+fi
+if [ ! -e sodShock.hdf5 ]
+then
+ echo "Generating initial conditions for the Traphic example..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../swift -R -t 4 traphic.yml 2>&1 | tee output.log
diff --git a/examples/Traphic/traphic.yml b/examples/Traphic/traphic.yml
new file mode 100644
index 0000000000000000000000000000000000000000..658d70bd82d4df87ea53394ce83412b2f6955dbb
--- /dev/null
+++ b/examples/Traphic/traphic.yml
@@ -0,0 +1,35 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1 # Grams
+ UnitLength_in_cgs: 1 # Centimeters
+ UnitVelocity_in_cgs: 1 # Centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 0.2 # The end time of the simulation (in internal units).
+ dt_min: 1e-7 # 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: traphic # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 0.2 # Time difference between consecutive outputs (in internal units)
+ compression: 1
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-2 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./traphic.hdf5 # The file to read
+
diff --git a/examples/main.c b/examples/main.c
index 773accc461fa26c819f0b603e7de689fa50e6cc6..b7647bbe3c073e7a2ac210a0d5f8d2a42cc9fea5 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -96,6 +96,7 @@ void print_help_message(void) {
"Set parameter value and overwrites values read from the parameters "
"file. Can be used more than once.");
printf(" %2s %14s %s\n", "-r", "", "Continue using restart files.");
+ printf(" %2s %14s %s\n", "-R", "", "Run with radiation.");
printf(" %2s %14s %s\n", "-s", "", "Run with hydrodynamics.");
printf(" %2s %14s %s\n", "-S", "", "Run with stars.");
printf(" %2s %14s %s\n", "-t", "{int}",
@@ -135,6 +136,7 @@ int main(int argc, char *argv[]) {
struct gpart *gparts = NULL;
struct gravity_props gravity_properties;
struct hydro_props hydro_properties;
+ struct rad_props radiation_properties;
struct part *parts = NULL;
struct phys_const prog_const;
struct sourceterms sourceterms;
@@ -196,6 +198,7 @@ int main(int argc, char *argv[]) {
int with_drift_all = 0;
int with_mpole_reconstruction = 0;
int with_structure_finding = 0;
+ int with_radiation = 0;
int verbose = 0;
int nr_threads = 1;
int with_verbose_timers = 0;
@@ -208,7 +211,7 @@ int main(int argc, char *argv[]) {
/* Parse the parameters */
int c;
- while ((c = getopt(argc, argv, "acCdDef:FgGhMn:o:P:rsSt:Tv:xy:Y:")) != -1)
+ while ((c = getopt(argc, argv, "acCdDef:FgGhMn:o:P:rRsSt:Tv:xy:Y:")) != -1)
switch (c) {
case 'a':
#if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA)
@@ -281,6 +284,9 @@ int main(int argc, char *argv[]) {
case 'r':
restart = 1;
break;
+ case 'R':
+ with_radiation = 1;
+ break;
case 's':
with_hydro = 1;
break;
@@ -369,9 +375,10 @@ int main(int argc, char *argv[]) {
if (myrank == 0) print_help_message();
return 1;
}
- if (!with_self_gravity && !with_hydro && !with_external_gravity) {
+ if (!with_self_gravity && !with_hydro && !with_external_gravity &&
+ !with_radiation) {
if (myrank == 0)
- printf("Error: At least one of -s, -g or -G must be chosen.\n");
+ printf("Error: At least one of -s, -R, -g or -G must be chosen.\n");
if (myrank == 0) print_help_message();
return 1;
}
@@ -459,6 +466,7 @@ int main(int argc, char *argv[]) {
message("sizeof(xpart) is %4zi bytes.", sizeof(struct xpart));
message("sizeof(spart) is %4zi bytes.", sizeof(struct spart));
message("sizeof(gpart) is %4zi bytes.", sizeof(struct gpart));
+ message("sizeof(rpart) is %4zi bytes.", sizeof(struct rpart));
message("sizeof(multipole) is %4zi bytes.", sizeof(struct multipole));
message("sizeof(grav_tensor) is %4zi bytes.", sizeof(struct grav_tensor));
message("sizeof(task) is %4zi bytes.", sizeof(struct task));
@@ -659,6 +667,9 @@ int main(int argc, char *argv[]) {
if (with_self_gravity)
gravity_props_init(&gravity_properties, params, &cosmo, with_cosmology);
+ /* Initialise the radiation properties */
+ if (with_radiation) radiation_props_init(&radiation_properties, params);
+
/* Read particles and space information from (GADGET) ICs */
char ICfileName[200] = "";
parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
@@ -706,11 +717,11 @@ int main(int argc, char *argv[]) {
dry_run);
#endif
#else
- read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
- &Ngpart, &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
- (with_external_gravity || with_self_gravity), with_stars,
- cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads,
- dry_run);
+ read_ic_single(
+ ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas, &Ngpart, &Nspart,
+ &periodic, &flag_entropy_ICs, with_hydro | with_radiation,
+ (with_external_gravity || with_self_gravity), with_stars, cleanup_h,
+ cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads, dry_run);
#endif
#endif
if (myrank == 0) {
@@ -784,7 +795,7 @@ int main(int argc, char *argv[]) {
if (myrank == 0) clocks_gettime(&tic);
space_init(&s, params, &cosmo, dim, parts, gparts, sparts, Ngas, Ngpart,
Nspart, periodic, replicate, generate_gas_in_ics,
- with_self_gravity, talking, dry_run);
+ with_self_gravity, with_radiation, talking, dry_run);
if (myrank == 0) {
clocks_gettime(&toc);
@@ -871,13 +882,14 @@ int main(int argc, char *argv[]) {
if (with_stars) engine_policies |= engine_policy_stars;
if (with_structure_finding)
engine_policies |= engine_policy_structure_finding;
+ if (with_radiation) engine_policies |= engine_policy_radiation;
/* Initialize the engine with the space and policies. */
if (myrank == 0) clocks_gettime(&tic);
engine_init(&e, &s, params, N_total[0], N_total[1], N_total[2],
engine_policies, talking, &reparttype, &us, &prog_const, &cosmo,
- &hydro_properties, &gravity_properties, &mesh, &potential,
- &cooling_func, &chemistry, &sourceterms);
+ &hydro_properties, &gravity_properties, &radiation_properties,
+ &mesh, &potential, &cooling_func, &chemistry, &sourceterms);
engine_config(0, &e, params, nr_nodes, myrank, nr_threads, with_aff,
talking, restart_file);
diff --git a/src/Makefile.am b/src/Makefile.am
index 0f61fb108d8d8acbd420c994266f4803a3f69d3e..e6834f4c5060e2534f33818201d8b6d49de5769c 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -48,7 +48,8 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \
gravity_softened_derivatives.h vector_power.h collectgroup.h hydro_space.h sort_part.h \
chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h space_getsid.h utilities.h \
- mesh_gravity.h cbrt.h velociraptor_interface.h swift_velociraptor_part.h outputlist.h
+ mesh_gravity.h cbrt.h velociraptor_interface.h swift_velociraptor_part.h outputlist.h \
+ radiation_properties.h radiation.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
@@ -61,7 +62,8 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
part_type.c xmf.c gravity_properties.c gravity.c \
collectgroup.c hydro_space.c equation_of_state.c \
chemistry.c cosmology.c restart.c mesh_gravity.c velociraptor_interface.c \
- outputlist.c
+ outputlist.c \
+ radiation_properties.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
@@ -152,7 +154,9 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
chemistry/EAGLE/chemistry.h \
chemistry/EAGLE/chemistry_io.h \
chemistry/EAGLE/chemistry_struct.h\
- chemistry/EAGLE/chemistry_iact.h
+ chemistry/EAGLE/chemistry_iact.h \
+ radiation/Traphic/rad_part.h \
+ radiation/Traphic/radiation.h
# Sources and flags for regular library
diff --git a/src/cell.c b/src/cell.c
index 05a1990dd5355877cadcd0526252a6266dddcf6b..b930fe17bd5b6312d946384ff70bcf92c010cd07 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -2380,6 +2380,28 @@ void cell_set_super_hydro(struct cell *c, struct cell *super_hydro) {
cell_set_super_hydro(c->progeny[k], super_hydro);
}
+/**
+ * @brief Set the super-cell pointers for all cells in a hierarchy.
+ *
+ * @param c The top-level #cell to play with.
+ * @param super_rad Pointer to the deepest cell with tasks in this part of the
+ * tree.
+ */
+void cell_set_super_radiation(struct cell *c, struct cell *super_rad) {
+
+ /* Are we in a cell with some kind of self/pair task ? */
+ if (super_rad == NULL && c->rad_density != NULL) super_rad = c;
+
+ /* Set the super-cell */
+ c->super_rad = super_rad;
+
+ /* Recurse */
+ if (c->split)
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL)
+ cell_set_super_radiation(c->progeny[k], super_rad);
+}
+
/**
* @brief Set the super-cell pointers for all cells in a hierarchy.
*
@@ -2424,6 +2446,8 @@ void cell_set_super_mapper(void *map_data, int num_elements, void *extra_data) {
(e->policy & engine_policy_external_gravity))
cell_set_super_gravity(c, NULL);
+ if (e->policy & engine_policy_radiation) cell_set_super_radiation(c, NULL);
+
/* Super-pointer for common operations */
cell_set_super(c, NULL);
}
diff --git a/src/cell.h b/src/cell.h
index 31d525b02b49463563b21bf5aa904a8b4301f989..096f44082c63a4898add8712c50c5bcf4329dc61 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -184,6 +184,9 @@ struct cell {
/*! Pointer to the #spart data. */
struct spart *sparts;
+ /*! Pointer to the #rpart data. */
+ struct rpart *rparts;
+
/*! Pointer for the sorted indices. */
struct entry *sort[13];
@@ -204,6 +207,10 @@ struct cell {
* tasks */
struct cell *super_gravity;
+ /*! Super cell, i.e. the highest-level parent clel that has radiation
+ * density pair/self tasks */
+ struct cell *super_rad;
+
/*! Linked list of the tasks computing this cell's hydro density. */
struct link *density;
@@ -216,6 +223,12 @@ struct cell {
/*! Linked list of the tasks computing this cell's gravity forces. */
struct link *grav;
+ /*! Linked list of the tasks computing this cell's radiation density. */
+ struct link *rad_density;
+
+ /*! Linked list of the tasks computing this cell's radiation transmission. */
+ struct link *rad_transmission;
+
/*! The task computing this cell's sorts. */
struct task *sorts;
diff --git a/src/cosmology.c b/src/cosmology.c
index d304354f3e5068137096c971fed693104924cdb3..65842176cdbbc068c5f247153e5585d6a9ff414e 100644
--- a/src/cosmology.c
+++ b/src/cosmology.c
@@ -374,7 +374,6 @@ void cosmology_init_tables(struct cosmology *c) {
GSL_INTEG_GAUSS61, space, &result, &abserr);
c->universe_age_at_present_day = result;
-
/* Update the times */
c->time_begin = cosmology_get_time_since_big_bang(c, c->a_begin);
c->time_end = cosmology_get_time_since_big_bang(c, c->a_end);
@@ -383,8 +382,7 @@ void cosmology_init_tables(struct cosmology *c) {
* Inverse t(a)
*/
- const double delta_t =
- (c->time_end - c->time_begin) / cosmology_table_length;
+ const double delta_t = (c->time_end - c->time_begin) / cosmology_table_length;
/* index in the time_interp_table */
int i_a = 0;
@@ -411,7 +409,7 @@ void cosmology_init_tables(struct cosmology *c) {
/* Compute interpolated scale factor */
double log_a = c->log_a_begin + scale * (c->log_a_end - c->log_a_begin) /
- cosmology_table_length;
+ cosmology_table_length;
c->scale_factor_interp_table[i_time] = exp(log_a) - c->a_begin;
}
diff --git a/src/engine.c b/src/engine.c
index d0668fa3c33f2a17ae112216d3b33d031376277e..026304edd691157c2ffcbbd129d6b1c891ae0c61 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -109,7 +109,8 @@ const char *engine_policy_names[] = {"none",
"cooling",
"sourceterms",
"stars",
- "structure finding"};
+ "structure finding",
+ "radiation"};
/** The rank of the engine as a global variable (for messages). */
int engine_rank;
@@ -297,6 +298,57 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c) {
}
}
+/**
+ * @brief Generate the radiation hierarchical tasks for a hierarchy of cells -
+ * i.e. all the O(Npart) tasks -- hydro version
+ *
+ * Tasks are only created here. The dependencies will be added later on.
+ *
+ * Note that there is no need to recurse below the super-cell. Note also
+ * that we only add tasks if the relevant particles are present in the cell.
+ *
+ * @param e The #engine.
+ * @param c The #cell.
+ */
+void engine_make_hierarchical_tasks_radiation(struct engine *e,
+ struct cell *c) {
+
+ struct scheduler *s = &e->sched;
+
+ /* Are we in a super-cell ? */
+ if (c->super_hydro == c) {
+
+ /* Add the sort task. */
+ c->sorts =
+ scheduler_addtask(s, task_type_sort, task_subtype_none, 0, 0, c, NULL);
+
+ /* Local tasks only... */
+ if (c->nodeID == e->nodeID) {
+
+ /* Generate the ghost tasks. */
+ c->ghost_in =
+ scheduler_addtask(s, task_type_ghost_in, task_subtype_none, 0,
+ /* implicit = */ 1, c, NULL);
+ c->ghost_out =
+ scheduler_addtask(s, task_type_ghost_out, task_subtype_none, 0,
+ /* implicit = */ 1, c, NULL);
+ engine_add_ghosts(e, c, c->ghost_in, c->ghost_out);
+
+ /* Generate the extra ghost task. */
+ c->extra_ghost = scheduler_addtask(s, task_type_extra_ghost,
+ task_subtype_none, 0, 0, c, NULL);
+ }
+
+ } else { /* We are above the super-cell so need to go deeper */
+
+ /* Recurse. */
+ if (c->split)
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL)
+ engine_make_hierarchical_tasks_radiation(e, c->progeny[k]);
+ }
+}
+
/**
* @brief Generate the hydro hierarchical tasks for a hierarchy of cells -
* i.e. all the O(Npart) tasks -- gravity version
@@ -396,6 +448,7 @@ void engine_make_hierarchical_tasks_mapper(void *map_data, int num_elements,
const int is_with_self_gravity = (e->policy & engine_policy_self_gravity);
const int is_with_external_gravity =
(e->policy & engine_policy_external_gravity);
+ const int is_with_radiation = (e->policy & engine_policy_radiation);
for (int ind = 0; ind < num_elements; ind++) {
struct cell *c = &((struct cell *)map_data)[ind];
@@ -406,6 +459,7 @@ void engine_make_hierarchical_tasks_mapper(void *map_data, int num_elements,
/* And the gravity stuff */
if (is_with_self_gravity || is_with_external_gravity)
engine_make_hierarchical_tasks_gravity(e, c);
+ if (is_with_radiation) engine_make_hierarchical_tasks_radiation(e, c);
}
}
@@ -2557,6 +2611,72 @@ void engine_make_external_gravity_tasks(struct engine *e) {
}
}
+/**
+ * @brief Constructs the top-level tasks for the radiation.
+ *
+ * @param e The #engine.
+ */
+void engine_make_first_radiation_tasks(struct engine *e) {
+
+ struct space *s = e->s;
+ struct scheduler *sched = &e->sched;
+ const int nodeID = e->nodeID;
+ const int *cdim = s->cdim;
+ struct cell *cells = s->cells_top;
+ const int nr_cells = s->nr_cells;
+
+ /* Loop through the elements, which are just byte offsets from NULL. */
+ for (int cid = 0; cid < nr_cells; ++cid) {
+
+ /* Get the cell index. */
+ const int i = cid / (cdim[1] * cdim[2]);
+ const int j = (cid / cdim[2]) % cdim[1];
+ const int k = cid % cdim[2];
+
+ /* Get the cell */
+ struct cell *ci = &cells[cid];
+
+ /* Skip cells without hydro particles */
+ if (ci->count == 0) continue;
+
+ /* If the cells is local build a self-interaction */
+ if (ci->nodeID == nodeID)
+ scheduler_addtask(sched, task_type_self, task_subtype_rad_density, 0, 0,
+ ci, NULL);
+
+ /* Now loop over all the neighbours of this cell */
+ for (int ii = -1; ii < 2; ii++) {
+ int iii = i + ii;
+ if (!s->periodic && (iii < 0 || iii >= cdim[0])) continue;
+ iii = (iii + cdim[0]) % cdim[0];
+ for (int jj = -1; jj < 2; jj++) {
+ int jjj = j + jj;
+ if (!s->periodic && (jjj < 0 || jjj >= cdim[1])) continue;
+ jjj = (jjj + cdim[1]) % cdim[1];
+ for (int kk = -1; kk < 2; kk++) {
+ int kkk = k + kk;
+ if (!s->periodic && (kkk < 0 || kkk >= cdim[2])) continue;
+ kkk = (kkk + cdim[2]) % cdim[2];
+
+ /* Get the neighbouring cell */
+ const int cjd = cell_getid(cdim, iii, jjj, kkk);
+ struct cell *cj = &cells[cjd];
+
+ /* Is that neighbour local and does it have particles ? */
+ if (cid >= cjd || cj->count == 0 ||
+ (ci->nodeID != nodeID && cj->nodeID != nodeID))
+ continue;
+
+ /* Construct the pair task */
+ const int sid = sortlistID[(kk + 1) + 3 * ((jj + 1) + 3 * (ii + 1))];
+ scheduler_addtask(sched, task_type_pair, task_subtype_rad_density,
+ sid, 0, ci, cj);
+ }
+ }
+ }
+ }
+}
+
/**
* @brief Constructs the top-level pair tasks for the first hydro loop over
* neighbours
@@ -2673,6 +2793,8 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
engine_addlink(e, &ci->grav, t);
} else if (t_subtype == task_subtype_external_grav) {
engine_addlink(e, &ci->grav, t);
+ } else if (t_subtype == task_subtype_rad_density) {
+ engine_addlink(e, &ci->rad_density, t);
}
/* Link pair tasks to cells. */
@@ -2686,6 +2808,9 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
} else if (t_subtype == task_subtype_grav) {
engine_addlink(e, &ci->grav, t);
engine_addlink(e, &cj->grav, t);
+ } else if (t_subtype == task_subtype_rad_density) {
+ engine_addlink(e, &ci->rad_density, t);
+ engine_addlink(e, &cj->rad_density, t);
}
#ifdef SWIFT_DEBUG_CHECKS
else if (t_subtype == task_subtype_external_grav) {
@@ -2703,6 +2828,8 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
engine_addlink(e, &ci->grav, t);
} else if (t_subtype == task_subtype_external_grav) {
engine_addlink(e, &ci->grav, t);
+ } else if (t_subtype == task_subtype_rad_density) {
+ engine_addlink(e, &ci->rad_density, t);
}
/* Link sub-pair tasks to cells. */
@@ -2716,6 +2843,9 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
} else if (t_subtype == task_subtype_grav) {
engine_addlink(e, &ci->grav, t);
engine_addlink(e, &cj->grav, t);
+ } else if (t_subtype == task_subtype_rad_density) {
+ engine_addlink(e, &ci->rad_density, t);
+ engine_addlink(e, &cj->rad_density, t);
}
#ifdef SWIFT_DEBUG_CHECKS
else if (t_subtype == task_subtype_external_grav) {
@@ -3174,6 +3304,162 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
}
}
+/**
+ * @brief Creates the dependency network for the radiation tasks of a given
+ * cell.
+ *
+ * @param sched The #scheduler.
+ * @param density The density task to link.
+ * @param force The force task to link.
+ * @param c The cell.
+ */
+static inline void engine_make_radiation_loops_dependencies(
+ struct scheduler *sched, struct task *rad_density,
+ struct task *rad_transmission, struct cell *c) {
+ /* density loop --> ghost --> transmission loop */
+ scheduler_addunlock(sched, rad_density, c->super_hydro->ghost_in);
+ scheduler_addunlock(sched, c->super_hydro->ghost_out, rad_transmission);
+}
+
+/**
+ * @brief Duplicates the first radiation loop and construct all the
+ * dependencies for the radiation parts
+ *
+ * This is done by looping over all the previously constructed tasks
+ * and adding another task involving the same cells but this time
+ * corresponding to the second radiation loop over neighbours.
+ * With all the relevant tasks for a given cell available, we construct
+ * all the dependencies for that cell.
+ */
+void engine_make_extra_radiation_tasks_mapper(void *map_data, int num_elements,
+ void *extra_data) {
+
+ struct engine *e = (struct engine *)extra_data;
+ struct scheduler *sched = &e->sched;
+ const int nodeID = e->nodeID;
+
+ for (int ind = 0; ind < num_elements; ind++) {
+ struct task *t = &((struct task *)map_data)[ind];
+
+ /* Self-interaction? */
+ if (t->type == task_type_self && t->subtype == task_subtype_rad_density) {
+
+ /* The self radiation density task has no dependencies */
+
+ /* Start by constructing the task for the second radiation loop */
+ struct task *t2 =
+ scheduler_addtask(sched, task_type_self,
+ task_subtype_rad_transmission, 0, 0, t->ci, NULL);
+
+ /* Add the link between the new loop and the cell */
+ engine_addlink(e, &t->ci->rad_transmission, t2);
+
+ /* Now, build all the dependencies for the hydro */
+ engine_make_radiation_loops_dependencies(sched, t, t2, t->ci);
+ /// if it doesn't work, you might want to look here
+ scheduler_addunlock(sched, t2, t->ci->super->end_force);
+ }
+
+ /* Otherwise, pair interaction? */
+ else if (t->type == task_type_pair &&
+ t->subtype == task_subtype_rad_density) {
+
+ /* Make all density tasks depend on the drift and the sorts. */
+ if (t->ci->nodeID == engine_rank)
+ scheduler_addunlock(sched, t->ci->super_hydro->drift_part, t);
+ scheduler_addunlock(sched, t->ci->super_hydro->sorts, t);
+ if (t->ci->super_hydro != t->cj->super_hydro) {
+ if (t->cj->nodeID == engine_rank)
+ scheduler_addunlock(sched, t->cj->super_hydro->drift_part, t);
+ scheduler_addunlock(sched, t->cj->super_hydro->sorts, t);
+ }
+
+ /* Start by constructing the task for the second hydro loop */
+ struct task *t2 = scheduler_addtask(
+ sched, task_type_pair, task_subtype_force, 0, 0, t->ci, t->cj);
+
+ /* Add the link between the new loop and both cells */
+ engine_addlink(e, &t->ci->force, t2);
+ engine_addlink(e, &t->cj->force, t2);
+
+ /* Now, build all the dependencies for the hydro for the cells */
+ /* that are local and are not descendant of the same super_hydro-cells */
+ if (t->ci->nodeID == nodeID) {
+ engine_make_radiation_loops_dependencies(sched, t, t2, t->ci);
+ scheduler_addunlock(sched, t2, t->ci->super->end_force);
+ }
+ if (t->cj->nodeID == nodeID) {
+ if (t->ci->super_hydro != t->cj->super_hydro)
+ engine_make_radiation_loops_dependencies(sched, t, t2, t->cj);
+ if (t->ci->super != t->cj->super)
+ scheduler_addunlock(sched, t2, t->cj->super->end_force);
+ }
+
+ }
+
+ /* Otherwise, sub-self interaction? */
+ else if (t->type == task_type_sub_self &&
+ t->subtype == task_subtype_density) {
+
+ /* Make all density tasks depend on the drift and sorts. */
+ scheduler_addunlock(sched, t->ci->super_hydro->drift_part, t);
+ scheduler_addunlock(sched, t->ci->super_hydro->sorts, t);
+
+ /* Start by constructing the task for the second hydro loop */
+ struct task *t2 =
+ scheduler_addtask(sched, task_type_sub_self, task_subtype_force,
+ t->flags, 0, t->ci, t->cj);
+
+ /* Add the link between the new loop and the cell */
+ engine_addlink(e, &t->ci->force, t2);
+
+ /* Now, build all the dependencies for the hydro for the cells */
+ /* that are local and are not descendant of the same super_hydro-cells */
+ if (t->ci->nodeID == nodeID) {
+ engine_make_radiation_loops_dependencies(sched, t, t2, t->ci);
+ scheduler_addunlock(sched, t2, t->ci->super->end_force);
+ }
+ }
+
+ /* Otherwise, sub-pair interaction? */
+ else if (t->type == task_type_sub_pair &&
+ t->subtype == task_subtype_density) {
+
+ /* Make all density tasks depend on the drift. */
+ if (t->ci->nodeID == engine_rank)
+ scheduler_addunlock(sched, t->ci->super_hydro->drift_part, t);
+ scheduler_addunlock(sched, t->ci->super_hydro->sorts, t);
+ if (t->ci->super_hydro != t->cj->super_hydro) {
+ if (t->cj->nodeID == engine_rank)
+ scheduler_addunlock(sched, t->cj->super_hydro->drift_part, t);
+ scheduler_addunlock(sched, t->cj->super_hydro->sorts, t);
+ }
+
+ /* Start by constructing the task for the second hydro loop */
+ struct task *t2 =
+ scheduler_addtask(sched, task_type_sub_pair, task_subtype_force,
+ t->flags, 0, t->ci, t->cj);
+
+ /* Add the link between the new loop and both cells */
+ engine_addlink(e, &t->ci->force, t2);
+ engine_addlink(e, &t->cj->force, t2);
+
+ /* Now, build all the dependencies for the hydro for the cells */
+ /* that are local and are not descendant of the same super_hydro-cells */
+ if (t->ci->nodeID == nodeID) {
+ engine_make_radiation_loops_dependencies(sched, t, t2, t->ci);
+ scheduler_addunlock(sched, t2, t->ci->super->end_force);
+ }
+ if (t->cj->nodeID == nodeID) {
+ if (t->ci->super_hydro != t->cj->super_hydro)
+ engine_make_radiation_loops_dependencies(sched, t, t2, t->cj);
+ if (t->ci->super != t->cj->super)
+ scheduler_addunlock(sched, t2, t->cj->super->end_force);
+ }
+ }
+ }
+}
+
/**
* @brief Fill the #space's task list.
*
@@ -3214,6 +3500,8 @@ void engine_maketasks(struct engine *e) {
if (e->policy & engine_policy_external_gravity)
engine_make_external_gravity_tasks(e);
+ if (e->policy & engine_policy_radiation) engine_make_first_radiation_tasks(e);
+
if (e->sched.nr_tasks == 0 && (s->nr_gparts > 0 || s->nr_parts > 0))
error("We have particles but no hydro or gravity tasks were created.");
@@ -3233,6 +3521,7 @@ void engine_maketasks(struct engine *e) {
#endif
const size_t self_grav_tasks_per_cell = 125;
const size_t ext_grav_tasks_per_cell = 1;
+ const size_t rad_tasks_per_cell = 27 * 2;
if (e->policy & engine_policy_hydro)
e->size_links += s->tot_cells * hydro_tasks_per_cell;
@@ -3240,6 +3529,8 @@ void engine_maketasks(struct engine *e) {
e->size_links += s->tot_cells * ext_grav_tasks_per_cell;
if (e->policy & engine_policy_self_gravity)
e->size_links += s->tot_cells * self_grav_tasks_per_cell;
+ if (e->policy & engine_policy_radiation)
+ e->size_links += s->tot_cells * rad_tasks_per_cell;
/* Allocate the new link list */
if ((e->links = (struct link *)malloc(sizeof(struct link) * e->size_links)) ==
@@ -3300,6 +3591,10 @@ void engine_maketasks(struct engine *e) {
threadpool_map(&e->threadpool, engine_make_extra_hydroloop_tasks_mapper,
sched->tasks, sched->nr_tasks, sizeof(struct task), 0, e);
+ if (e->policy & engine_policy_radiation)
+ threadpool_map(&e->threadpool, engine_make_extra_radiation_tasks_mapper,
+ sched->tasks, sched->nr_tasks, sizeof(struct task), 0, e);
+
if (e->verbose)
message("Making extra hydroloop tasks took %.3f %s.",
clocks_from_ticks(getticks() - tic2), clocks_getunit());
@@ -3891,6 +4186,10 @@ int engine_estimate_nr_tasks(struct engine *e) {
if (e->policy & engine_policy_stars) {
n1 += 2;
}
+ if (e->policy & engine_policy_radiation) {
+ /* 2 x 13 x neighbours + 2 x self + 2 x ghost */
+ n1 += 30;
+ }
#ifdef WITH_MPI
@@ -5717,6 +6016,7 @@ void engine_unpin(void) {
* @param cosmo The #cosmology used for this run.
* @param hydro The #hydro_props used for this run.
* @param gravity The #gravity_props used for this run.
+ * @param radiation The #rad_props used for this run.
* @param mesh The #pm_mesh used for the long-range periodic forces.
* @param potential The properties of the external potential.
* @param cooling_func The properties of the cooling function.
@@ -5729,7 +6029,8 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
const struct unit_system *internal_units,
const struct phys_const *physical_constants,
struct cosmology *cosmo, const struct hydro_props *hydro,
- struct gravity_props *gravity, struct pm_mesh *mesh,
+ struct gravity_props *gravity,
+ const struct rad_props *radiation, struct pm_mesh *mesh,
const struct external_potential *potential,
const struct cooling_function_data *cooling_func,
const struct chemistry_global_data *chemistry,
@@ -5795,6 +6096,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
e->cosmology = cosmo;
e->hydro_properties = hydro;
e->gravity_properties = gravity;
+ e->radiation_properties = radiation;
e->mesh = mesh;
e->external_potential = potential;
e->cooling_func = cooling_func;
diff --git a/src/engine.h b/src/engine.h
index cfd656712bdea84484101cf7e83795f795200f5f..5cedbe0506f2d1df929f097bfd96387a4a4c61e4 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -71,9 +71,10 @@ enum engine_policy {
engine_policy_cooling = (1 << 13),
engine_policy_sourceterms = (1 << 14),
engine_policy_stars = (1 << 15),
- engine_policy_structure_finding = (1 << 16)
+ engine_policy_structure_finding = (1 << 16),
+ engine_policy_radiation = (1 << 17)
};
-#define engine_maxpolicy 16
+#define engine_maxpolicy 17
extern const char *engine_policy_names[];
/**
@@ -326,6 +327,9 @@ struct engine {
/* Properties of the self-gravity scheme */
struct gravity_props *gravity_properties;
+ /* Properties of the radiation scheme */
+ const struct rad_props *radiation_properties;
+
/* The mesh used for long-range gravity forces */
struct pm_mesh *mesh;
@@ -390,7 +394,8 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
const struct unit_system *internal_units,
const struct phys_const *physical_constants,
struct cosmology *cosmo, const struct hydro_props *hydro,
- struct gravity_props *gravity, struct pm_mesh *mesh,
+ struct gravity_props *gravity,
+ const struct rad_props *radiation, struct pm_mesh *mesh,
const struct external_potential *potential,
const struct cooling_function_data *cooling_func,
const struct chemistry_global_data *chemistry,
diff --git a/src/part.h b/src/part.h
index 145bf2111771d8ad254affb213b93b7ac829f1a6..74a3bef11d2887a7261e25463aadafad7b292e9d 100644
--- a/src/part.h
+++ b/src/part.h
@@ -40,6 +40,7 @@
#define xpart_align 128
#define spart_align 128
#define gpart_align 128
+#define rpart_align 128
/* Import the right hydro particle definition */
#if defined(MINIMAL_SPH)
@@ -88,6 +89,9 @@
/* Import the right star particle definition */
#include "./stars/Default/star_part.h"
+/* Import the right radiation particle definition */
+#include "./radiation/Traphic/rad_part.h"
+
void part_relink_gparts_to_parts(struct part *parts, size_t N,
ptrdiff_t offset);
void part_relink_gparts_to_sparts(struct spart *sparts, size_t N,
diff --git a/src/radiation.h b/src/radiation.h
new file mode 100644
index 0000000000000000000000000000000000000000..451227eeef02ad3c3538db3391cba04ead79c59d
--- /dev/null
+++ b/src/radiation.h
@@ -0,0 +1,32 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_RADIATION_H
+#define SWIFT_RADIATION_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local headers. */
+#include "const.h"
+#include "inline.h"
+#include "part.h"
+
+#include "./radiation/Traphic/radiation.h"
+
+#endif
diff --git a/src/radiation/Traphic/rad_part.h b/src/radiation/Traphic/rad_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..f824fb3ea3031e269871b5a206abed1deb8f2bdc
--- /dev/null
+++ b/src/radiation/Traphic/rad_part.h
@@ -0,0 +1,47 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_TRAPHIC_RAD_PART_H
+#define SWIFT_TRAPHIC_RAD_PART_H
+
+/* Some standard headers. */
+#include <stdlib.h>
+
+/**
+ * @brief Particle fields for the radiation particles.
+ */
+struct rpart {
+
+ /*! Particle position. */
+ double x[3];
+
+ /*! Particle smoothing length. */
+ float h;
+
+ /*! Particle density. */
+ double density;
+
+ /*! Particle hydrogen neutral fraction. */
+ double hydrogen_neutral_fraction;
+
+ /*! Particle received ionising luminosity. */
+ double ionising_luminosity;
+
+} SWIFT_STRUCT_ALIGN;
+
+#endif /* SWIFT_TRAPHIC_RAD_PART_H */
diff --git a/src/radiation/Traphic/radiation.h b/src/radiation/Traphic/radiation.h
new file mode 100644
index 0000000000000000000000000000000000000000..db4bb14df5b7a1021b06f2baa10a8409edf663b2
--- /dev/null
+++ b/src/radiation/Traphic/radiation.h
@@ -0,0 +1,41 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_TRAPHIC_RADIATION_H
+#define SWIFT_TRAPHIC_RADIATION_H
+
+/**
+ * @brief Initialise the rpart based on the original gas particle.
+ *
+ * @param p Original gas particle.
+ * @param rp #rpart to initialise.
+ */
+__attribute__((always_inline)) INLINE static void radiation_first_init_part(
+ struct part* p, struct rpart* rp) {
+
+ rp->x[0] = p->x[0];
+ rp->x[1] = p->x[1];
+ rp->x[2] = p->x[2];
+ rp->h = p->h;
+
+ rp->density = 1.;
+ rp->hydrogen_neutral_fraction = 1.;
+ rp->ionising_luminosity = 0.;
+}
+
+#endif /* SWIFT_TRAPHIC_RADIATION_H */
diff --git a/src/radiation_properties.c b/src/radiation_properties.c
new file mode 100644
index 0000000000000000000000000000000000000000..b6b243d2e9403b2f7f1ca64d755373d1e0037e16
--- /dev/null
+++ b/src/radiation_properties.c
@@ -0,0 +1,61 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* This object's header. */
+#include "radiation_properties.h"
+
+/* Standard headers */
+#include <float.h>
+#include <math.h>
+
+/* Local headers. */
+#include "common_io.h"
+#include "error.h"
+
+void radiation_props_init(struct rad_props *p, struct swift_params *params) {}
+
+void radiation_props_print(const struct rad_props *p) {}
+
+#if defined(HAVE_HDF5)
+void radiation_props_print_snapshot(hid_t h_grpgrav,
+ const struct rad_props *p) {}
+#endif
+
+/**
+ * @brief Write a rad_props struct to the given FILE as a stream of bytes.
+ *
+ * @param p the struct
+ * @param stream the file stream
+ */
+void radiation_props_struct_dump(const struct rad_props *p, FILE *stream) {
+ restart_write_blocks((void *)p, sizeof(struct rad_props), 1, stream,
+ "radiation", "rad props");
+}
+
+/**
+ * @brief Restore a rad_props struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param p the struct
+ * @param stream the file stream
+ */
+void radiation_props_struct_restore(struct rad_props *p, FILE *stream) {
+ restart_read_blocks((void *)p, sizeof(struct rad_props), 1, stream, NULL,
+ "rad props");
+}
diff --git a/src/radiation_properties.h b/src/radiation_properties.h
new file mode 100644
index 0000000000000000000000000000000000000000..376fbe556d6c2c4cf36a9c2732e6732994471c2f
--- /dev/null
+++ b/src/radiation_properties.h
@@ -0,0 +1,52 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_RADIATION_PROPERTIES
+#define SWIFT_RADIATION_PROPERTIES
+
+/* Config parameters. */
+#include "../config.h"
+
+#if defined(HAVE_HDF5)
+#include <hdf5.h>
+#endif
+
+/* Local includes. */
+#include "restart.h"
+
+/* Forward declarations */
+struct swift_params;
+
+/**
+ * @brief Contains all the constants and parameters of the radiation transfer
+ * scheme.
+ */
+struct rad_props {};
+
+void radiation_props_print(const struct rad_props *p);
+void radiation_props_init(struct rad_props *p, struct swift_params *params);
+
+#if defined(HAVE_HDF5)
+void radiation_props_print_snapshot(hid_t h_grpsph, const struct rad_props *p);
+#endif
+
+/* Dump/restore. */
+void radiation_props_struct_dump(const struct rad_props *p, FILE *stream);
+void radiation_props_struct_restore(struct rad_props *p, FILE *stream);
+
+#endif /* SWIFT_RADIATION_PROPERTIES */
diff --git a/src/scheduler.c b/src/scheduler.c
index 4974884651b02db57d851493a4fc8fe343483a05..0cf898b3a812cd3a858095fbfea1e48f9d714006 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -985,6 +985,8 @@ void scheduler_splittasks_mapper(void *map_data, int num_elements,
scheduler_splittask_gravity(t, s);
} else if (t->type == task_type_grav_mesh) {
/* For future use */
+ } else if (t->subtype == task_subtype_rad_density) {
+ scheduler_splittask_hydro(t, s);
} else {
#ifdef SWIFT_DEBUG_CHECKS
error("Unexpected task sub-type");
diff --git a/src/space.c b/src/space.c
index 6f98e788e9625c1cc872f59c58a8bf87b7b2cfa8..541bcb70e314c7d3de8b7035855f233882570a4d 100644
--- a/src/space.c
+++ b/src/space.c
@@ -54,6 +54,7 @@
#include "memswap.h"
#include "minmax.h"
#include "multipole.h"
+#include "radiation.h"
#include "restart.h"
#include "sort_part.h"
#include "stars.h"
@@ -2446,6 +2447,14 @@ void space_first_init_parts_mapper(void *restrict map_data, int count,
p[k].ti_kick = 0;
#endif
}
+
+ /* Initialize the rparts */
+ if (s->radiation) {
+ struct rpart *restrict rp = s->rparts + delta;
+ for (int k = 0; k < count; k++) {
+ radiation_first_init_part(&p[k], &rp[k]);
+ }
+ }
}
/**
@@ -2688,6 +2697,7 @@ void space_convert_quantities(struct space *s, int verbose) {
* @param replicate How many replications along each direction do we want?
* @param generate_gas_in_ics Are we generating gas particles from the gparts?
* @param self_gravity flag whether we are doing gravity or not?
+ * @param radiation Flag whether we are doing radiation or not?
* @param verbose Print messages to stdout or not.
* @param dry_run If 1, just initialise stuff, don't do anything with the parts.
*
@@ -2701,7 +2711,7 @@ void space_init(struct space *s, struct swift_params *params,
struct part *parts, struct gpart *gparts, struct spart *sparts,
size_t Npart, size_t Ngpart, size_t Nspart, int periodic,
int replicate, int generate_gas_in_ics, int self_gravity,
- int verbose, int dry_run) {
+ int radiation, int verbose, int dry_run) {
/* Clean-up everything */
bzero(s, sizeof(struct space));
@@ -2712,6 +2722,7 @@ void space_init(struct space *s, struct swift_params *params,
s->dim[2] = dim[2];
s->periodic = periodic;
s->gravity = self_gravity;
+ s->radiation = radiation;
s->nr_parts = Npart;
s->size_parts = Npart;
s->parts = parts;
@@ -2890,6 +2901,14 @@ void space_init(struct space *s, struct swift_params *params,
bzero(s->xparts, Npart * sizeof(struct xpart));
}
+ /* Allocate the extra radiation parts array for the gas particles. */
+ if (radiation) {
+ if (posix_memalign((void **)&s->rparts, rpart_align,
+ Npart * sizeof(struct rpart)) != 0)
+ error("Failed to allocate rparts.");
+ bzero(s->rparts, Npart * sizeof(struct rpart));
+ }
+
hydro_space_init(&s->hs, s);
/* Init the space lock. */
diff --git a/src/space.h b/src/space.h
index e3173ece1e2749a3afb8072b179150587a100a82..2bba805c5fa5788a260c14e914dd6deca2479c1c 100644
--- a/src/space.h
+++ b/src/space.h
@@ -82,6 +82,9 @@ struct space {
/*! Are we doing gravity? */
int gravity;
+ /*! Are we doing radiation? */
+ int radiation;
+
/*! Width of the top-level cells. */
double width[3];
@@ -145,6 +148,9 @@ struct space {
/*! The s-particle data (cells have pointers to this). */
struct spart *sparts;
+ /*! The r-particle data (cells have pointers to this). */
+ struct rpart *rparts;
+
/*! The top-level FFT task */
struct task *grav_top_level;
@@ -207,7 +213,7 @@ void space_init(struct space *s, struct swift_params *params,
struct part *parts, struct gpart *gparts, struct spart *sparts,
size_t Npart, size_t Ngpart, size_t Nspart, int periodic,
int replicate, int generate_gas_in_ics, int gravity,
- int verbose, int dry_run);
+ int radiation, int verbose, int dry_run);
void space_sanitize(struct space *s);
void space_map_cells_pre(struct space *s, int full,
void (*fun)(struct cell *c, void *data), void *data);
diff --git a/src/swift.h b/src/swift.h
index e10938addb99956c202b3e4dd2b0592b580fa948..df6b36816f56206d4b15ac85bf8eafdeef2e6084 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -59,6 +59,8 @@
#include "potential.h"
#include "profiler.h"
#include "queue.h"
+#include "radiation.h"
+#include "radiation_properties.h"
#include "restart.h"
#include "runner.h"
#include "scheduler.h"
diff --git a/src/task.h b/src/task.h
index 072d3979ce04990aaef46c5cc5eb0b8c62fdc860..5e412b9efd60cd0b96acc8b723517abb7fb67620 100644
--- a/src/task.h
+++ b/src/task.h
@@ -85,6 +85,8 @@ enum task_subtypes {
task_subtype_gpart,
task_subtype_multipole,
task_subtype_spart,
+ task_subtype_rad_density,
+ task_subtype_rad_transmission,
task_subtype_count
} __attribute__((packed));