diff --git a/examples/DiscPatch/HydroStatic/dynamic.pro b/examples/DiscPatch/HydroStatic/dynamic.pro
index c02c65fe418e84cdd62978dbddcf5a641fa4c156..00ee3f7a8d2dc435be2093af959efd2c49903637 100644
--- a/examples/DiscPatch/HydroStatic/dynamic.pro
+++ b/examples/DiscPatch/HydroStatic/dynamic.pro
@@ -8,7 +8,8 @@ iplot = 1 ; if iplot = 1, make plot of E/Lz conservation, else, simply compare f
@physunits
indir = './'
-basefile = 'Disk-Patch-dynamic_'
+;basefile = 'Disc-Patch-dynamic_'
+basefile = 'Disc-Patch_'
; set properties of potential
uL = phys.pc ; unit of length
@@ -16,18 +17,27 @@ uM = phys.msun ; unit of mass
uV = 1d5 ; unit of velocity
; properties of patch
-surface_density = 10.
+surface_density = 100. ; surface density of all mass, which generates the gravitational potential
scale_height = 100.
-z_disk = 200.;
+z_disk = 200. ;
+fgas = 0.1 ; gas fraction
gamma = 5./3.
; derived units
constG = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
pcentre = [0.,0.,z_disk] * pc / uL
utherm = !pi * constG * surface_density * scale_height / (gamma-1.)
+temp = (utherm*uV^2)*phys.m_h/phys.kb
soundspeed = sqrt(gamma * (gamma-1.) * utherm)
t_dyn = sqrt(scale_height / (constG * surface_density))
-
+rho0 = fgas*(surface_density)/(2.*scale_height)
+print,' dynamical time = ',t_dyn,' = ',t_dyn*UL/uV/(1d6*phys.yr),' Myr'
+print,' thermal energy per unit mass = ',utherm
+print,' central density = ',rho0,' = ',rho0*uM/uL^3/m_h,' particles/cm^3'
+print,' central temperature = ',temp
+lambda = 2 * !pi * phys.G^1.5 * (scale_height*uL)^1.5 * (surface_density * uM/uL^2)^0.5 * phys.m_h^2 / (gamma-1) / fgas
+print,' lambda = ',lambda
+stop
;
infile = indir + basefile + '*'
spawn,'ls -1 '+infile,res
diff --git a/examples/DiscPatch/HydroStatic/makeIC.py b/examples/DiscPatch/HydroStatic/makeIC.py
index 48cc578658a9520477d40bf504e3eb7c3c8e5164..e2846d82a8cfa8bf08de83632b19ae2e7818f3c1 100644
--- a/examples/DiscPatch/HydroStatic/makeIC.py
+++ b/examples/DiscPatch/HydroStatic/makeIC.py
@@ -56,9 +56,10 @@ print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs
# parameters of potential
-surface_density = 10.
+surface_density = 100. # surface density of all mass, which generates the gravitational potential
scale_height = 100.
gamma = 5./3.
+fgas = 0.1 # gas fraction
# derived units
const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
@@ -131,7 +132,7 @@ h = glass_h[0:numGas]
numGas = numpy.shape(pos)[0]
# compute furthe properties of ICs
-column_density = surface_density * numpy.tanh(boxSize/2./scale_height)
+column_density = fgas * surface_density * numpy.tanh(boxSize/2./scale_height)
enclosed_mass = column_density * boxSize * boxSize
pmass = enclosed_mass / numGas
meanrho = enclosed_mass / boxSize**3
diff --git a/examples/Feedback/feedback.pro b/examples/Feedback/feedback.pro
new file mode 100644
index 0000000000000000000000000000000000000000..02d616fc82f0aeb7011d022d13db9d1d1030e89c
--- /dev/null
+++ b/examples/Feedback/feedback.pro
@@ -0,0 +1,24 @@
+base = 'Feedback'
+inf = 'Feedback_005.hdf5'
+
+blast = [5.650488e-01, 5.004371e-01, 5.010494e-01] ; location of blast
+pos = h5rd(inf,'PartType0/Coordinates')
+vel = h5rd(inf,'PartType0/Velocities')
+rho = h5rd(inf,'PartType0/Density')
+utherm = h5rd(inf,'PartType0/InternalEnergy')
+
+; shift to centre
+for ic=0,2 do pos[ic,*] = pos[ic,*] - blast[ic]
+
+;; distance from centre
+dist = fltarr(n_elements(rho))
+for ic=0,2 do dist = dist + pos[ic,*]^2
+dist = sqrt(dist)
+
+; radial velocity
+vr = fltarr(n_elements(rho))
+for ic=0,2 do vr = vr + pos[ic,*]*vel[ic,*]
+vr = vr / dist
+
+;
+end
diff --git a/examples/Feedback/feedback.yml b/examples/Feedback/feedback.yml
new file mode 100644
index 0000000000000000000000000000000000000000..9ff7eea325b43fe3c670d16bdbd4ccc66f33333e
--- /dev/null
+++ b/examples/Feedback/feedback.yml
@@ -0,0 +1,44 @@
+# 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: 5e-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-4 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: Feedback # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 1e-2 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-3 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ max_smoothing_length: 0.1 # Maximal smoothing length allowed (in internal units).
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./Feedback.hdf5 # The file to read
+
+# Parameters for feedback
+
+SN:
+ time: 0.001 # time the SN explodes (internal units)
+ energy: 1.0 # energy of the explosion (internal units)
+ x: 0.5 # x-position of explostion (internal units)
+ y: 0.5 # y-position of explostion (internal units)
+ z: 0.5 # z-position of explostion (internal units)
diff --git a/examples/Feedback/makeIC.py b/examples/Feedback/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..bd1081a9c275616038f5fa4e3eb943c36cb4c3eb
--- /dev/null
+++ b/examples/Feedback/makeIC.py
@@ -0,0 +1,109 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ #
+ # 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
+import sys
+from numpy import *
+
+# Generates a swift IC file containing a cartesian distribution of particles
+# at a constant density and pressure in a cubic box
+
+# Parameters
+periodic= 1 # 1 For periodic box
+boxSize = 1.
+L = int(sys.argv[1]) # Number of particles along one axis
+rho = 1. # Density
+P = 1.e-6 # Pressure
+gamma = 5./3. # Gas adiabatic index
+eta = 1.2349 # 48 ngbs with cubic spline kernel
+fileName = "Feedback.hdf5"
+
+#---------------------------------------------------
+numPart = L**3
+mass = boxSize**3 * rho / numPart
+internalEnergy = P / ((gamma - 1.)*rho)
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+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
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+#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")
+
+v = zeros((numPart, 3))
+ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+ds[()] = v
+v = zeros(1)
+
+m = full((numPart, 1), mass)
+ds = grp.create_dataset('Masses', (numPart,1), 'f')
+ds[()] = m
+m = zeros(1)
+
+h = full((numPart, 1), eta * boxSize / L)
+ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
+ds[()] = h
+h = zeros(1)
+
+u = full((numPart, 1), internalEnergy)
+ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
+ds[()] = u
+u = zeros(1)
+
+
+ids = linspace(0, numPart, numPart, endpoint=False).reshape((numPart,1))
+ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
+ds[()] = ids + 1
+x = ids % L;
+y = ((ids - x) / L) % L;
+z = (ids - x - L * y) / L**2;
+coords = zeros((numPart, 3))
+coords[:,0] = z[:,0] * boxSize / L + boxSize / (2*L)
+coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
+coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
+ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = coords
+
+file.close()
diff --git a/examples/main.c b/examples/main.c
index 895bb9003e04a4723b13ea1eaa360c5aa44e7228..a8abb938398dbcd7fa22578cbae29ed58e8c8ffe 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -76,6 +76,7 @@ void print_help_message() {
"Overwrite the CPU frequency (Hz) to be used for time measurements");
printf(" %2s %8s %s\n", "-g", "",
"Run with an external gravitational potential");
+ printf(" %2s %8s %s\n", "-F", "", "Run with feedback ");
printf(" %2s %8s %s\n", "-G", "", "Run with self-gravity");
printf(" %2s %8s %s\n", "-n", "{int}",
"Execute a fixed number of time steps. When unset use the time_end "
@@ -147,6 +148,7 @@ int main(int argc, char *argv[]) {
int nsteps = -2;
int with_cosmology = 0;
int with_external_gravity = 0;
+ int with_sourceterms = 0;
int with_cooling = 0;
int with_self_gravity = 0;
int with_hydro = 0;
@@ -159,7 +161,7 @@ int main(int argc, char *argv[]) {
/* Parse the parameters */
int c;
- while ((c = getopt(argc, argv, "acCdDef:gGhn:st:v:y:")) != -1) switch (c) {
+ while ((c = getopt(argc, argv, "acCdDef:FgGhn:st:v:y:")) != -1) switch (c) {
case 'a':
with_aff = 1;
break;
@@ -185,6 +187,9 @@ int main(int argc, char *argv[]) {
return 1;
}
break;
+ case 'F':
+ with_sourceterms = 1;
+ break;
case 'g':
with_external_gravity = 1;
break;
@@ -443,6 +448,11 @@ int main(int argc, char *argv[]) {
if (with_cooling) cooling_init(params, &us, &prog_const, &cooling_func);
if (with_cooling && myrank == 0) cooling_print(&cooling_func);
+ /* Initialise the feedback properties */
+ struct sourceterms sourceterms;
+ if (with_sourceterms) sourceterms_init(params, &us, &sourceterms);
+ if (with_sourceterms && myrank == 0) sourceterms_print(&sourceterms);
+
/* Construct the engine policy */
int engine_policies = ENGINE_POLICY | engine_policy_steal;
if (with_drift_all) engine_policies |= engine_policy_drift_all;
@@ -451,13 +461,14 @@ int main(int argc, char *argv[]) {
if (with_external_gravity) engine_policies |= engine_policy_external_gravity;
if (with_cosmology) engine_policies |= engine_policy_cosmology;
if (with_cooling) engine_policies |= engine_policy_cooling;
+ if (with_sourceterms) engine_policies |= engine_policy_sourceterms;
/* Initialize the engine with the space and policies. */
if (myrank == 0) clocks_gettime(&tic);
struct engine e;
engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, with_aff,
engine_policies, talking, &us, &prog_const, &hydro_properties,
- &potential, &cooling_func);
+ &potential, &cooling_func, &sourceterms);
if (myrank == 0) {
clocks_gettime(&toc);
message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
diff --git a/src/Makefile.am b/src/Makefile.am
index 32c6f0a27200c8c416543a291241fb2a336e9037..5859c59d9ddbc7e701145c40f8f0a35b92a075c3 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -43,7 +43,8 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
engine.h swift.h serial_io.h timers.h debug.h scheduler.h proxy.h parallel_io.h \
common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h parser.h \
physical_constants.h physical_constants_cgs.h potential.h version.h \
- hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h
+ hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h \
+ sourceterms_struct.h
# Common source files
@@ -52,7 +53,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
units.c common_io.c single_io.c multipole.c version.c map.c \
kernel_hydro.c tools.c part.c partition.c clocks.c parser.c \
physical_constants.c potential.c hydro_properties.c \
- runner_doiact_fft.c threadpool.c cooling.c
+ runner_doiact_fft.c threadpool.c cooling.c sourceterms.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
@@ -62,6 +63,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h
gravity.h gravity_io.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
+ sourceterms.h \
hydro.h hydro_io.h \
hydro/Minimal/hydro.h hydro/Minimal/hydro_iact.h hydro/Minimal/hydro_io.h \
hydro/Minimal/hydro_debug.h hydro/Minimal/hydro_part.h \
diff --git a/src/cell.h b/src/cell.h
index 8c64ebe3e6e82e279284a56b20f25d487ff0b836..576d2db695720b3aa1b49b89128496a7ca12a6bb 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -163,6 +163,9 @@ struct cell {
/* Task for cooling */
struct task *cooling;
+ /* Task for source terms */
+ struct task *sourceterms;
+
/* Number of tasks that are associated with this cell. */
int nr_tasks;
diff --git a/src/const.h b/src/const.h
index afad73211546cf81cbb0e6a1179569e588952fec..ef26525b7865bac4c0164f7678398ca307fc1d67 100644
--- a/src/const.h
+++ b/src/const.h
@@ -96,6 +96,10 @@
//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
//#define EXTERNAL_POTENTIAL_DISC_PATCH
+/* Source terms */
+#define SOURCETERMS_NONE
+//#define SOURCETERMS_SN_FEEDBACK
+
/* Cooling properties */
#define COOLING_NONE
//#define COOLING_CONST_DU
diff --git a/src/engine.c b/src/engine.c
index c1c37dde79b7b3891f4e4c743184bd7a0c92c386..1ebc83972e747cefdbdd84cf2a23b3effbd53985 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -68,7 +68,7 @@
#include "units.h"
#include "version.h"
-const char *engine_policy_names[15] = {"none",
+const char *engine_policy_names[16] = {"none",
"rand",
"steal",
"keep",
@@ -82,7 +82,8 @@ const char *engine_policy_names[15] = {"none",
"external_gravity",
"cosmology_integration",
"drift_all",
- "cooling"};
+ "cooling",
+ "sourceterms"};
/** The rank of the engine as a global variable (for messages). */
int engine_rank;
@@ -128,6 +129,9 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
const int is_fixdt = (e->policy & engine_policy_fixdt);
const int is_hydro = (e->policy & engine_policy_hydro);
const int is_with_cooling = (e->policy & engine_policy_cooling);
+ const int is_with_sourceterms =
+ (e->policy & engine_policy_sourceterms) == engine_policy_sourceterms;
+
/* Are we in a super-cell ? */
if (c->super == c) {
@@ -164,6 +168,10 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
if (is_with_cooling)
c->cooling = scheduler_addtask(s, task_type_cooling, task_subtype_none,
0, 0, c, NULL, 0);
+ /* add source terms */
+ if (is_with_sourceterms)
+ c->sourceterms = scheduler_addtask(s, task_type_sourceterms,
+ task_subtype_none, 0, 0, c, NULL, 0);
}
} else { /* We are above the super-cell so need to go deeper */
@@ -1783,12 +1791,18 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
}
#endif
}
-
- /* Cooling tasks should depend on kick and does not unlock anything since
- it is the last task*/
+ /* Cooling tasks should depend on kick and unlock sourceterms */
else if (t->type == task_type_cooling) {
scheduler_addunlock(sched, t->ci->kick, t);
}
+ /* source terms depend on cooling if performed, else on kick. It is the last
+ task */
+ else if (t->type == task_type_sourceterms) {
+ if (e->policy == engine_policy_cooling)
+ scheduler_addunlock(sched, t->ci->cooling, t);
+ else
+ scheduler_addunlock(sched, t->ci->kick, t);
+ }
}
}
@@ -2824,6 +2838,11 @@ void engine_step(struct engine *e) {
mask |= 1 << task_type_cooling;
}
+ /* Add the tasks corresponding to sourceterms to the masks */
+ if (e->policy & engine_policy_sourceterms) {
+ mask |= 1 << task_type_sourceterms;
+ }
+
/* Add MPI tasks if need be */
if (e->policy & engine_policy_mpi) {
@@ -3167,6 +3186,7 @@ void engine_unpin() {
* @param hydro The #hydro_props used for this run.
* @param potential The properties of the external potential.
* @param cooling_func The properties of the cooling function.
+ * @param sourceterms The properties of the source terms function.
*/
void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
@@ -3175,7 +3195,8 @@ void engine_init(struct engine *e, struct space *s,
const struct phys_const *physical_constants,
const struct hydro_props *hydro,
const struct external_potential *potential,
- const struct cooling_function_data *cooling_func) {
+ const struct cooling_function_data *cooling_func,
+ struct sourceterms *sourceterms) {
/* Clean-up everything */
bzero(e, sizeof(struct engine));
@@ -3228,6 +3249,7 @@ void engine_init(struct engine *e, struct space *s,
e->hydro_properties = hydro;
e->external_potential = potential;
e->cooling_func = cooling_func;
+ e->sourceterms = sourceterms;
e->parameter_file = params;
engine_rank = nodeID;
diff --git a/src/engine.h b/src/engine.h
index d36914af611723bc4d496d5c2dc68050eea6ffe6..ec9e16553c3172f22e9bfe60971163488947a47e 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -44,6 +44,7 @@
#include "potential.h"
#include "runner.h"
#include "scheduler.h"
+#include "sourceterms_struct.h"
#include "space.h"
#include "task.h"
#include "units.h"
@@ -65,6 +66,7 @@ enum engine_policy {
engine_policy_cosmology = (1 << 11),
engine_policy_drift_all = (1 << 12),
engine_policy_cooling = (1 << 13),
+ engine_policy_sourceterms = (1 << 14)
};
extern const char *engine_policy_names[];
@@ -207,6 +209,9 @@ struct engine {
/* Properties of the cooling scheme */
const struct cooling_function_data *cooling_func;
+ /* Properties of source terms */
+ struct sourceterms *sourceterms;
+
/* The (parsed) parameter file */
const struct swift_params *parameter_file;
};
@@ -223,7 +228,8 @@ void engine_init(struct engine *e, struct space *s,
const struct phys_const *physical_constants,
const struct hydro_props *hydro,
const struct external_potential *potential,
- const struct cooling_function_data *cooling);
+ const struct cooling_function_data *cooling,
+ struct sourceterms *sourceterms);
void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
unsigned int submask);
void engine_prepare(struct engine *e, int nodrift);
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 282f7f481944b8637831f565ee73ed1184a66e76..0bdf56fe844309998db1fad4cf9581edb6b88305 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -37,6 +37,7 @@
#include "equation_of_state.h"
#include "hydro_properties.h"
#include "kernel_hydro.h"
+#include "minmax.h"
/**
* @brief Returns the internal energy of a particle
@@ -116,8 +117,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
* @brief Modifies the thermal state of a particle to the imposed internal
* energy
*
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
+ * updated.
*
* @param p The particle
* @param u The new internal energy
@@ -133,17 +135,24 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
/* Compute the new sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+ /* Update the signal velocity */
+ const float v_sig_old = p->force.v_sig;
+ const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+ const float v_sig = max(v_sig_old, v_sig_new);
+
const float rho_inv = 1.f / p->rho;
p->force.soundspeed = soundspeed;
p->force.P_over_rho2 = pressure * rho_inv * rho_inv;
+ p->force.v_sig = v_sig;
}
/**
* @brief Modifies the thermal state of a particle to the imposed entropy
*
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
+ * updated.
*
* @param p The particle
* @param S The new entropy
@@ -159,10 +168,16 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
/* Compute the new sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+ /* Update the signal velocity */
+ const float v_sig_old = p->force.v_sig;
+ const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+ const float v_sig = max(v_sig_old, v_sig_new);
+
const float rho_inv = 1.f / p->rho;
p->force.soundspeed = soundspeed;
p->force.P_over_rho2 = pressure * rho_inv * rho_inv;
+ p->force.v_sig = v_sig;
}
/**
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 42dc2aa3d8899ce01d1b7f0bc55f5faa43e45a68..3015f26c6bad7f006e8cda16695c750ff40d74df 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -39,6 +39,7 @@
#include "equation_of_state.h"
#include "hydro_properties.h"
#include "kernel_hydro.h"
+#include "minmax.h"
/**
* @brief Returns the internal energy of a particle
@@ -124,8 +125,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
* @brief Modifies the thermal state of a particle to the imposed internal
* energy
*
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Internal energy, pressure, sound-speed and signal velocity
+ * will be updated.
*
* @param p The particle
* @param u The new internal energy
@@ -141,15 +143,22 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
/* Compute the new sound speed */
const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
+ /* Update the signal velocity */
+ const float v_sig_old = p->force.v_sig;
+ const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+ const float v_sig = max(v_sig_old, v_sig_new);
+
p->force.soundspeed = soundspeed;
p->force.pressure = pressure;
+ p->force.v_sig = v_sig;
}
/**
* @brief Modifies the thermal state of a particle to the imposed entropy
*
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Internal energy, pressure, sound-speed and signal velocity
+ * will be updated.
*
* @param p The particle
* @param S The new entropy
@@ -165,8 +174,14 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
/* Compute the new sound speed */
const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
+ /* Update the signal velocity */
+ const float v_sig_old = p->force.v_sig;
+ const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+ const float v_sig = max(v_sig_old, v_sig_new);
+
p->force.soundspeed = soundspeed;
p->force.pressure = pressure;
+ p->force.v_sig = v_sig;
}
/**
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index a5e0d41804f5b03cc5e06dd452a022ba5eedb73f..ede16fe5181a546921db7acf8ffcbc130c5e91c5 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -37,6 +37,7 @@
#include "equation_of_state.h"
#include "hydro_properties.h"
#include "kernel_hydro.h"
+#include "minmax.h"
/**
* @brief Returns the internal energy of a particle
@@ -116,8 +117,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
* @brief Modifies the thermal state of a particle to the imposed internal
* energy
*
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
+ * updated.
*
* @param p The particle
* @param u The new internal energy
@@ -133,17 +135,25 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
/* Compute the sound speed from the pressure*/
const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
+
+ /* Update the signal velocity */
+ const float v_sig_old = p->force.v_sig;
+ const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+ const float v_sig = max(v_sig_old, v_sig_new);
+
const float rho_bar_inv = 1.f / p->rho_bar;
p->force.soundspeed = soundspeed;
p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
+ p->force.v_sig = v_sig;
}
/**
* @brief Modifies the thermal state of a particle to the imposed entropy
*
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
+ * updated.
*
* @param p The particle
* @param S The new entropy
@@ -159,10 +169,17 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
/* Compute the sound speed from the pressure*/
const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
+
+ /* Update the signal velocity */
+ const float v_sig_old = p->force.v_sig;
+ const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+ const float v_sig = max(v_sig_old, v_sig_new);
+
const float rho_bar_inv = 1.f / p->rho_bar;
p->force.soundspeed = soundspeed;
p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
+ p->force.v_sig = v_sig;
}
/**
diff --git a/src/runner.c b/src/runner.c
index fe81947a49194398682f621b0acefaa6ad249572..06cef10cde74f034dd9b15fe7ef124f267341a25 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -53,6 +53,7 @@
#include "kick.h"
#include "minmax.h"
#include "scheduler.h"
+#include "sourceterms.h"
#include "space.h"
#include "task.h"
#include "timers.h"
@@ -111,6 +112,42 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
#include "runner_doiact_fft.h"
#include "runner_doiact_grav.h"
+/**
+ * @brief Perform source terms
+ *
+ * @param r runner task
+ * @param c cell
+ * @param timer 1 if the time is to be recorded.
+ */
+void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) {
+ const int count = c->count;
+ const double cell_min[3] = {c->loc[0], c->loc[1], c->loc[2]};
+ const double cell_width[3] = {c->width[0], c->width[1], c->width[2]};
+ struct sourceterms *sourceterms = r->e->sourceterms;
+ const int dimen = 3;
+
+ TIMER_TIC;
+
+ /* Recurse? */
+ if (c->split) {
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) runner_do_sourceterms(r, c->progeny[k], 0);
+ return;
+ }
+
+ if (count > 0) {
+
+ /* do sourceterms in this cell? */
+ const int incell =
+ sourceterms_test_cell(cell_min, cell_width, sourceterms, dimen);
+ if (incell == 1) {
+ sourceterms_apply(r, sourceterms, c);
+ }
+ }
+
+ if (timer) TIMER_TOC(timer_dosource);
+}
+
/**
* @brief Calculate gravity acceleration from external potential
*
@@ -1345,6 +1382,9 @@ void *runner_main(void *data) {
case task_type_cooling:
runner_do_cooling(r, t->ci, 1);
break;
+ case task_type_sourceterms:
+ runner_do_sourceterms(r, t->ci, 1);
+ break;
default:
error("Unknown task type.");
}
diff --git a/src/sourceterms.c b/src/sourceterms.c
new file mode 100644
index 0000000000000000000000000000000000000000..2a5f326688b52b0e4abfde8dca823d765e0e8a5e
--- /dev/null
+++ b/src/sourceterms.c
@@ -0,0 +1,60 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local includes. */
+#include "const.h"
+#include "hydro.h"
+#include "parser.h"
+#include "units.h"
+
+/* This object's header. */
+#include "sourceterms.h"
+
+/**
+ * @brief Initialises the sourceterms
+ *
+ * @param parameter_file The parsed parameter file
+ * @param us The current internal system of units
+ * @param source the structure that has all the source term properties
+ */
+void sourceterms_init(const struct swift_params* parameter_file,
+ struct UnitSystem* us, struct sourceterms* source) {
+#ifdef SOURCETERMS_SN_FEEDBACK
+ supernova_init(parameter_file, us, source);
+#endif /* SOURCETERMS_SN_FEEDBACK */
+};
+
+/**
+ * @brief Prints the properties of the source terms to stdout
+ * @param source the structure that has all the source term properties
+ */
+void sourceterms_print(struct sourceterms* source) {
+#ifdef SOURCETERMS_NONE
+ error(" no sourceterms defined yet you ran with -F");
+#ifdef SOURCETERMS_SN_FEEDBACK
+#error "can't have sourceterms when defined SOURCETERMS_NONE"
+#endif
+#endif
+#ifdef SOURCETERMS_SN_FEEDBACK
+ supernova_print(source);
+#endif /* SOURCETERMS_SN_FEEDBACK */
+};
diff --git a/src/sourceterms.h b/src/sourceterms.h
new file mode 100644
index 0000000000000000000000000000000000000000..b6fb6d18094b6cd42717eb6931b3dd5f8947433c
--- /dev/null
+++ b/src/sourceterms.h
@@ -0,0 +1,74 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * 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_SOURCETERMS_H
+#define SWIFT_SOURCETERMS_H
+
+/**
+ * @file src/sourceterms.h
+ * @brief Branches between the different sourceterms functions.
+ */
+
+#include "./const.h"
+#include "runner.h"
+
+#ifdef SOURCETERMS_SN_FEEDBACK
+#include "sourceterms/sn_feedback/sn_feedback_struct.h"
+#endif
+
+/* So far only one model here */
+struct sourceterms {
+#ifdef SOURCETERMS_SN_FEEDBACK
+ struct supernova_struct supernova;
+#endif
+};
+#ifdef SOURCETERMS_SN_FEEDBACK
+#include "sourceterms/sn_feedback/sn_feedback.h"
+#endif
+
+void sourceterms_init(const struct swift_params* parameter_file,
+ struct UnitSystem* us, struct sourceterms* source);
+void sourceterms_print(struct sourceterms* source);
+/**
+ * @file src/sourceterm.h
+ * @brief Routines related to source terms
+ * @param cell_min: corner of cell to test
+ * @param cell_width: width of cell to test
+ * @param sourceterms: properties of source terms to test
+ * @param dimen: dimensionality of the problem
+ *
+ * This routine tests whether a source term should be applied to this cell
+ * return: 1 if yes, return: 0 if no
+ */
+
+__attribute__((always_inline)) INLINE static int sourceterms_test_cell(
+ const double cell_min[], const double cell_width[],
+ struct sourceterms* sourceterms, const int dimen) {
+#ifdef SOURCETERMS_SN_FEEDBACK
+ return supernova_feedback_test_cell(cell_min, cell_width, sourceterms, dimen);
+#endif
+ return 0;
+};
+
+__attribute__((always_inline)) INLINE static void sourceterms_apply(
+ struct runner* r, struct sourceterms* sourceterms, struct cell* c) {
+#ifdef SOURCETERMS_SN_FEEDBACK
+ supernova_feedback_apply(r, sourceterms, c);
+#endif
+};
+#endif /* SWIFT_SOURCETERMS_H */
diff --git a/src/sourceterms/sn_feedback/sn_feedback.h b/src/sourceterms/sn_feedback/sn_feedback.h
new file mode 100644
index 0000000000000000000000000000000000000000..667c79cb1eef42f7fde79c43059d1baa5c61268e
--- /dev/null
+++ b/src/sourceterms/sn_feedback/sn_feedback.h
@@ -0,0 +1,192 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ *
+ * 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_SN_FEEDBACK_H
+#define SWIFT_SN_FEEDBACK_H
+#include <float.h>
+/* Config parameters. */
+#include "../config.h"
+
+#include "engine.h"
+#include "equation_of_state.h"
+#include "hydro.h"
+#include "runner.h"
+#include "timestep.h"
+
+/**
+ * @file src/sourceterms/sn_feedback.h
+ *
+ * @brief Routines related to sourceterms (supernova feedback): determine if
+ * feedback occurs in this cell
+ *
+ * @param cell_min: corner of cell to test
+ * @param cell_width: width of cell to test
+ * @param sourceterms: properties of source terms to test
+ * @param dimen: dimensionality of the problem
+ *
+ * This routine tests whether a source term should be applied to this cell
+ * return: 1 if yes, return: 0 if no
+ */
+__attribute__((always_inline)) INLINE static int supernova_feedback_test_cell(
+ const double cell_min[], const double cell_width[],
+ struct sourceterms* sourceterms, const int dimen) {
+ if (sourceterms->supernova.status == supernova_is_done) return 0;
+
+ const double location[3] = {sourceterms->supernova.x,
+ sourceterms->supernova.y,
+ sourceterms->supernova.z};
+ for (int i = 0; i < dimen; i++) {
+ if (cell_min[i] > location[i]) return 0;
+ if ((cell_min[i] + cell_width[i]) <= location[i]) return 0;
+ };
+ return 1;
+};
+
+/**
+ * @file src/sourceterms/sn_feedback.h
+ *
+ * @brief Routines related to source terms (supernova feedback): perform
+ * feedback in this cell
+ * @param r: the runner
+ * @param sourceterms the structure describing the source terms properties
+ * @param c the cell to apply feedback to
+ *
+ * This routine heats an individual particle (p), increasing its thermal energy
+ * per unit mass
+ * by supernova energy / particle mass.
+ */
+__attribute__((always_inline)) INLINE static void supernova_feedback_apply(
+ struct runner* restrict r, struct sourceterms* restrict sourceterms,
+ struct cell* restrict c) {
+
+ const int count = c->count;
+ struct part* restrict parts = c->parts;
+ struct xpart* restrict xparts = c->xparts;
+ const double timeBase = r->e->timeBase;
+ const int ti_current = r->e->ti_current;
+
+ /* inject SN energy into the particle with highest id in this cell if it is
+ * active */
+ int imax = 0;
+ struct part* restrict p_sn = NULL;
+ struct xpart* restrict xp_sn = NULL;
+
+ for (int i = 0; i < count; i++) {
+
+ /* Get a direct pointer on the part. */
+ struct part* restrict p = &parts[i];
+ if (p->id > imax) {
+ imax = p->id;
+ p_sn = p;
+ xp_sn = &xparts[i];
+ }
+ }
+
+ /* Is this part within the time step? */
+ if (p_sn->ti_begin == ti_current) {
+
+ /* Does this time step straddle the feedback injection time? */
+ const float t_begin = p_sn->ti_begin * timeBase;
+ const float t_end = p_sn->ti_end * timeBase;
+ if (t_begin <= sourceterms->supernova.time &&
+ t_end > sourceterms->supernova.time) {
+
+ /* store old time step */
+ const int dti_old = p_sn->ti_end - p_sn->ti_begin;
+
+ /* add supernova feedback */
+ const float u_old = hydro_get_internal_energy(p_sn, 0);
+ const float ent_old = hydro_get_entropy(p_sn, 0.0);
+ const float u_new =
+ u_old + sourceterms->supernova.energy / hydro_get_mass(p_sn);
+ hydro_set_internal_energy(p_sn, u_new);
+ const float u_set = hydro_get_internal_energy(p_sn, 0.0);
+ const float ent_set = hydro_get_entropy(p_sn, 0.0);
+ message(
+ " applied super nova, time = %e, location= %e %e %e velocity= %e %e "
+ "%e",
+ ti_current * timeBase, p_sn->x[0], p_sn->x[1], p_sn->x[2], p_sn->v[0],
+ p_sn->v[1], p_sn->v[2]);
+ message(
+ " injected SN energy in particle = %lld, increased energy from %e to "
+ "%e and is notw %e, entropy from %e to %e",
+ p_sn->id, u_old, u_new, u_set, ent_old, ent_set);
+
+ /* label supernova as done */
+ sourceterms->supernova.status = supernova_is_done;
+
+ /* update timestep if new time step shorter than old time step */
+ const int dti = get_part_timestep(p_sn, xp_sn, r->e);
+ if (dti < dti_old) {
+ p_sn->ti_end = p_sn->ti_begin + dti;
+ message(" changed timestep from %d to %d", dti_old, dti);
+
+ /* apply simple time-step limiter on all particles in same cell:
+ */
+ int i_limit = 0;
+ for (int i = 0; i < count; i++) {
+ struct part* restrict p = &parts[i];
+ const int dti_old = p->ti_end - p->ti_begin;
+ if (dti_old > 2 * dti) {
+ i_limit++;
+ const int dti_new = 2 * dti;
+ p->ti_end = p->ti_begin + dti_new;
+ message(" old step = %d new step = %d", dti_old, dti_new);
+ } else
+ message(" old step = %d", dti_old);
+ }
+ message(" count= %d limited timestep of %d particles ", count, i_limit);
+ } /* end of limiter */
+ error("end");
+ }
+ }
+};
+
+/**
+ * @file src/sourceterms/sn_feedback.h
+ *
+ * @brief Routine to initialise supernova feedback
+ * @param parameterfile: the parse parmeter file
+ * @param us: the unit system in use
+ * @param sourceterms the structure describing the source terms properties
+ *
+ * This routine heats an individual particle (p), increasing its thermal energy
+ * per unit mass
+ * by supernova energy / particle mass.
+ */
+
+__attribute__((always_inline)) INLINE static void supernova_init(
+ const struct swift_params* parameter_file, struct UnitSystem* us,
+ struct sourceterms* source) {
+ source->supernova.time = parser_get_param_double(parameter_file, "SN:time");
+ source->supernova.energy =
+ parser_get_param_double(parameter_file, "SN:energy");
+ source->supernova.x = parser_get_param_double(parameter_file, "SN:x");
+ source->supernova.y = parser_get_param_double(parameter_file, "SN:y");
+ source->supernova.z = parser_get_param_double(parameter_file, "SN:z");
+ source->supernova.status = supernova_is_not_done;
+}
+__attribute__((always_inline)) INLINE static void supernova_print(
+ struct sourceterms* source) {
+ message(
+ " Single SNe of energy= %e will explode at time= %e at location "
+ "(%e,%e,%e)",
+ source->supernova.energy, source->supernova.time, source->supernova.x,
+ source->supernova.y, source->supernova.z);
+}
+#endif /* SWIFT_SN_FEEDBACK_H */
diff --git a/src/sourceterms/sn_feedback/sn_feedback_struct.h b/src/sourceterms/sn_feedback/sn_feedback_struct.h
new file mode 100644
index 0000000000000000000000000000000000000000..2b4dbe97b0a687cc48d0ba7e12b105f306cb7e96
--- /dev/null
+++ b/src/sourceterms/sn_feedback/sn_feedback_struct.h
@@ -0,0 +1,45 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+/**
+ * @file src/sourceterms/sn_feedback_struct.h
+ * @brief Routines related to source terms (feedback)
+ *
+ * enumeration type that sets if supernova explosion is done (is_done) or still
+ * needs doing (is_not_done)
+ */
+#ifndef SWIFT_SN_FEEDBACK_STRUCT_H
+#define SWIFT_SN_FEEDBACK_STRUCT_H
+enum supernova_status { supernova_is_done, supernova_is_not_done };
+
+/**
+ * @file src/sourceterms/sn_feedback_struct.h
+ * @brief Routines related to source terms (feedback)
+ *
+ * The structure that describes the source term (supernova feedback)
+ * It specifies the time, energy and location of the desired supernova
+ * explosion, and a status (supernova_is_done/supernova_is_not_done)
+ * that records the status of the supernova
+ */
+struct supernova_struct {
+ double time;
+ double energy;
+ double x, y, z;
+ enum supernova_status status;
+};
+#endif SWIFT_SN_FEEDBACK_STRUCT_H
diff --git a/src/sourceterms_struct.h b/src/sourceterms_struct.h
new file mode 100644
index 0000000000000000000000000000000000000000..b3c38986db52d72df825fda97b36c985dff922b6
--- /dev/null
+++ b/src/sourceterms_struct.h
@@ -0,0 +1,26 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * 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_SOURCETERMS_STRUCT_H
+#define SWIFT_SOURCETERMS_STRUCT_H
+#include "./const.h"
+#ifdef SOURCETERMS_SN_FEEDBACK
+#include "sourceterms/sn_feedback/sn_feedback_struct.h"
+#endif
+
+#endif /* SWIFT_SOURCETERMS_STRUCT_H */
diff --git a/src/swift.h b/src/swift.h
index 3119adfd4b95639a84137e243c6e971d80ef3825..58745f3111fc3f29490d4591ee549907392cb87e 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -49,6 +49,7 @@
#include "scheduler.h"
#include "serial_io.h"
#include "single_io.h"
+#include "sourceterms.h"
#include "space.h"
#include "task.h"
#include "timers.h"
diff --git a/src/task.c b/src/task.c
index 30e90f90a420a2efc742ec2c8d8454bf01d42e86..473df4aecb13823f4d49a0d2c7d64c528ab1c3c1 100644
--- a/src/task.c
+++ b/src/task.c
@@ -50,7 +50,8 @@
const char *taskID_names[task_type_count] = {
"none", "sort", "self", "pair", "sub_self", "sub_pair",
"init", "ghost", "extra_ghost", "kick", "kick_fixdt", "send",
- "recv", "grav_gather_m", "grav_fft", "grav_mm", "grav_up", "cooling"};
+ "recv", "grav_gather_m", "grav_fft", "grav_mm", "grav_up", "cooling",
+ "sourceterms"};
const char *subtaskID_names[task_subtype_count] = {
"none", "density", "gradient", "force", "grav", "external_grav", "tend"};
@@ -117,6 +118,7 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
case task_type_ghost:
case task_type_extra_ghost:
case task_type_cooling:
+ case task_type_sourceterms:
return task_action_part;
break;
diff --git a/src/task.h b/src/task.h
index efd3596144cc741bb84597cdfbecd402ba02bc33..f840c0b4b8e807dce28f6f13479dbdf4995ab66d 100644
--- a/src/task.h
+++ b/src/task.h
@@ -54,6 +54,7 @@ enum task_types {
task_type_grav_mm,
task_type_grav_up,
task_type_cooling,
+ task_type_sourceterms,
task_type_count
} __attribute__((packed));
diff --git a/src/timers.h b/src/timers.h
index b93a34df4d90251d4686afaef0be51b36dc9a25e..d75508afd88cf2fd57d06f9530bff8607d79d127 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -45,6 +45,7 @@ enum {
timer_dopair_grav_pm,
timer_dopair_grav_pp,
timer_dograv_external,
+ timer_dosource,
timer_dosub_self_density,
timer_dosub_self_gradient,
timer_dosub_self_force,
diff --git a/src/timestep.h b/src/timestep.h
index 599fb4762e11b08fc942fb02acbbf1970f477de4..db52911ec1e8fbf31f35e8877e0a7ae7ba5ee478 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -24,8 +24,8 @@
/* Local headers. */
#include "const.h"
+#include "cooling.h"
#include "debug.h"
-
/**
* @brief Compute a valid integer time-step form a given time-step
*