diff --git a/examples/Feedback/feedback.pro b/examples/Feedback/feedback.pro
index 93305120e00d54a377bf0a0ee75bcc2fe7586ab0..02d616fc82f0aeb7011d022d13db9d1d1030e89c 100644
--- a/examples/Feedback/feedback.pro
+++ b/examples/Feedback/feedback.pro
@@ -1,7 +1,7 @@
base = 'Feedback'
-inf = 'Feedback_003.hdf5'
+inf = 'Feedback_005.hdf5'
-blast = [0.5,0.5,0.5] ; location of blast
+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')
@@ -10,7 +10,7 @@ utherm = h5rd(inf,'PartType0/InternalEnergy')
; shift to centre
for ic=0,2 do pos[ic,*] = pos[ic,*] - blast[ic]
-;; distnace from centre
+;; distance from centre
dist = fltarr(n_elements(rho))
for ic=0,2 do dist = dist + pos[ic,*]^2
dist = sqrt(dist)
@@ -19,4 +19,6 @@ dist = sqrt(dist)
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
index 74d2445973db9be6a955efe9be5be86705c25bc1..9ff7eea325b43fe3c670d16bdbd4ccc66f33333e 100644
--- a/examples/Feedback/feedback.yml
+++ b/examples/Feedback/feedback.yml
@@ -37,8 +37,8 @@ InitialConditions:
# Parameters for feedback
SN:
- time: 0.001
- energy: 1.0
- x: 0.5
- y: 0.5
- z: 0.5
+ 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
index d494a99e456d11019aae1450233125ffaf07f7ad..bd1081a9c275616038f5fa4e3eb943c36cb4c3eb 100644
--- a/examples/Feedback/makeIC.py
+++ b/examples/Feedback/makeIC.py
@@ -2,6 +2,7 @@
# 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
diff --git a/src/engine.c b/src/engine.c
index f6ccd4a6b8ac3f223defa6a6439628ca5ea421cf..c320cbc68654ff6a3ab9742c94163f03a8712376 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1799,14 +1799,15 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
/* Cooling tasks should depend on kick and unlock sourceterms */
else if (t->type == task_type_cooling) {
scheduler_addunlock(sched, t->ci->kick, t);
- scheduler_addunlock(sched, t->ci->sourceterms, t);
- }
- /* source terms depend on cooling if performed, else on kick. It is the last task */
+ scheduler_addunlock(sched, t->ci->sourceterms, 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);
+ if (e->policy == engine_policy_cooling)
+ scheduler_addunlock(sched, t->ci->cooling, t);
+ else
+ scheduler_addunlock(sched, t->ci->kick, t);
}
}
}
diff --git a/src/runner.c b/src/runner.c
index b77e6978f0b4acbaa2cd6fe7ff3f4a2fc81a3400..558f29a46281c27ae4b24fab150699d20d698b09 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -133,7 +133,6 @@ void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) {
sourceterms->supernova.z};
const int dimen = 3;
const double timeBase = r->e->timeBase;
- const float CFL_condition = r->e->hydro_properties->CFL_condition;
TIMER_TIC;
@@ -145,95 +144,80 @@ void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) {
}
/* is supernova still active? */
- if(sourceterms->supernova.status == supernova_is_not_done)
- {
- /* does cell contain explosion? */
- if (count > 0) {
- const int incell = is_in_cell(cell_min, cell_width, location, dimen);
- if (incell == 1) {
-
- /* inject SN energy into particle with highest id, 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) {
-
- /* add supernova feedback */
- const int orig_dti = get_part_timestep(p_sn, xp_sn, r->e);
- const float ci_old = hydro_get_soundspeed(p_sn, 0.0); /* sound speed */
- const float u_old = hydro_get_internal_energy(p_sn, 0.0);
- do_supernova_feedback(sourceterms, p_sn);
- /* updating u does not propagate corerctly in Gizmo branch: needs updating /
-
- /* label supernova as done */
- sourceterms->supernova.status = supernova_is_done;
- message(" applied super nova, time = %d", ti_current);
- /* recompute time step */
- const float u = hydro_get_internal_energy(p_sn, 0.0);
- const float ci = sqrt(hydro_gamma_minus_one * hydro_gamma * u);
- /* following needs updating to be amde consisten
- const float ci = hydro_get_soundspeed((p_sn, 0.0);
- */
- const float max_timestep = CFL_condition * p_sn->h / ci; /* new maximum timestep */
- const int old_dti = p_sn->ti_end - p_sn->ti_begin;
- const float old_timestep = (p_sn->ti_end - p_sn->ti_begin) * timeBase; /* current time step */
- float sn_timestep = old_timestep;
- message(" old c=%e old u= %e new c= %e new u= %e", ci_old, u_old, ci, u);
- message(" h= %e cfl= %e", p_sn->h, CFL_condition);
- message(" old dt= %e new dt = %e", old_timestep, max_timestep);
- if(max_timestep < old_timestep)
- {
- sn_timestep = max_timestep;
- const int new_dti =
- get_integer_timestep(max_timestep, p_sn->ti_begin, p_sn->ti_end, r->e->timeBase_inv);
- p_sn->ti_end = p_sn->ti_begin + new_dti;
- /* apply new time-step */
- p_sn->ti_end = p_sn->ti_begin + new_dti;
- message(" changed timestep from %d %d to %d", orig_dti, old_dti, new_dti);
-
- /* apply simple time-step limiter on all particles in this cell: */
- int i_limit=0;
- for (int i=0; i<count; i++)
- {
- struct part *restrict p = &parts[i];
- const float old_dt = (p->ti_end - p->ti_begin) * timeBase;
- if(old_dt > 2*sn_timestep)
- {
- i_limit++;
- const int new_dti =
- get_integer_timestep(max_timestep, p->ti_begin, p->ti_end, r->e->timeBase_inv);
- p->ti_end = p->ti_begin + new_dti;
- message(" old step = %e new step = %e", old_dt, new_dti * r->e->timeBase);
- }
- }
- message(" limited timestep of %d particles ",i_limit);
-
- }
- error(" check! ");
- }
- }
- }
- }
- }
+ if (sourceterms->supernova.status == supernova_is_not_done) {
+ /* does cell contain explosion? */
+ if (count > 0) {
+ const int incell = is_in_cell(cell_min, cell_width, location, dimen);
+ if (incell == 1) {
+
+ /* inject SN energy into particle with highest id, 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 */
+ do_supernova_feedback(sourceterms, p_sn);
+
+ /* label supernova as done */
+ sourceterms->supernova.status = supernova_is_done;
+ message(" applied super nova, time = %d, location= %e %e %e",
+ ti_current, p_sn->x[0], p_sn->x[1], p_sn->x[2]);
+ message(" applied super nova, velocity = %e %e %e", p_sn->v[0],
+ p_sn->v[1], p_sn->v[2]);
+ error("end");
+
+ /* 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);
+ }
+ // error(" check! ");
+ }
+ }
+ }
+ }
+ }
if (timer) TIMER_TOC(timer_dosource);
}
diff --git a/src/sourceterms.c b/src/sourceterms.c
index 6d424b8390c0497bfadaa0370b350285e6c41a40..13dac2c3f65b00903f692540de8f547b4556658f 100644
--- a/src/sourceterms.c
+++ b/src/sourceterms.c
@@ -36,7 +36,7 @@
* @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) {
+ struct UnitSystem* us, struct sourceterms* source) {
#ifdef SN_FEEDBACK
source->supernova.time = parser_get_param_double(parameter_file, "SN:time");
diff --git a/src/sourceterms.h b/src/sourceterms.h
index 15e26cc708a4c26a161ccc5600971e9cccb0628f..cac53e52df1667ab967b26fb2b98bcde8157aeb1 100644
--- a/src/sourceterms.h
+++ b/src/sourceterms.h
@@ -31,7 +31,7 @@ struct sourceterms {
};
void sourceterms_init(const struct swift_params* parameter_file,
- struct UnitSystem* us, struct sourceterms* source);
+ struct UnitSystem* us, struct sourceterms* source);
void sourceterms_print(struct sourceterms* source);
#ifdef SN_FEEDBACK
#include "sourceterms/sn_feedback/sn_feedback.h"
diff --git a/src/sourceterms/sn_feedback/sn_feedback.h b/src/sourceterms/sn_feedback/sn_feedback.h
index 4dd413e77ef4571ff0a9b80de31a662817c88ecd..e58f687767a6c8fbecb6569fd6445edb6a5244b7 100644
--- a/src/sourceterms/sn_feedback/sn_feedback.h
+++ b/src/sourceterms/sn_feedback/sn_feedback.h
@@ -38,7 +38,8 @@ __attribute__((always_inline)) INLINE static int is_in_cell(
* @param sourceterms the structure describing the source terms properties
* @param p the particle to apply feedback to
*
- * This routine heats an individual particle (p), increasing its thermal energy per unit mass
+ * 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 do_supernova_feedback(
@@ -49,7 +50,9 @@ __attribute__((always_inline)) INLINE static void do_supernova_feedback(
hydro_set_internal_energy(p, u_new);
const float u_set = hydro_get_internal_energy(p, 0.0);
message(" unew = %e %e s= %e", u_new, u_set, p->entropy);
- message(" injected SN energy in particle = %lld, increased energy from %e to %e, check= %e", p->id, u_old, u_new, u_set);
-
+ message(
+ " injected SN energy in particle = %lld, increased energy from %e to %e, "
+ "check= %e",
+ p->id, u_old, u_new, u_set);
};
#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
index 37190236ddfea6c8afaf5cf75436c6bf533913e6..2b4dbe97b0a687cc48d0ba7e12b105f306cb7e96 100644
--- a/src/sourceterms/sn_feedback/sn_feedback_struct.h
+++ b/src/sourceterms/sn_feedback/sn_feedback_struct.h
@@ -20,21 +20,20 @@
* @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)
+ * 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
-};
+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)
+ * 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 {