Commit b401d0d2 authored by Tom Theuns's avatar Tom Theuns
Browse files

merged with master

parent 6b6defc1
base = 'Feedback' 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') pos = h5rd(inf,'PartType0/Coordinates')
vel = h5rd(inf,'PartType0/Velocities') vel = h5rd(inf,'PartType0/Velocities')
rho = h5rd(inf,'PartType0/Density') rho = h5rd(inf,'PartType0/Density')
...@@ -10,7 +10,7 @@ utherm = h5rd(inf,'PartType0/InternalEnergy') ...@@ -10,7 +10,7 @@ utherm = h5rd(inf,'PartType0/InternalEnergy')
; shift to centre ; shift to centre
for ic=0,2 do pos[ic,*] = pos[ic,*] - blast[ic] for ic=0,2 do pos[ic,*] = pos[ic,*] - blast[ic]
;; distnace from centre ;; distance from centre
dist = fltarr(n_elements(rho)) dist = fltarr(n_elements(rho))
for ic=0,2 do dist = dist + pos[ic,*]^2 for ic=0,2 do dist = dist + pos[ic,*]^2
dist = sqrt(dist) dist = sqrt(dist)
...@@ -19,4 +19,6 @@ dist = sqrt(dist) ...@@ -19,4 +19,6 @@ dist = sqrt(dist)
vr = fltarr(n_elements(rho)) vr = fltarr(n_elements(rho))
for ic=0,2 do vr = vr + pos[ic,*]*vel[ic,*] for ic=0,2 do vr = vr + pos[ic,*]*vel[ic,*]
vr = vr / dist vr = vr / dist
;
end end
...@@ -37,8 +37,8 @@ InitialConditions: ...@@ -37,8 +37,8 @@ InitialConditions:
# Parameters for feedback # Parameters for feedback
SN: SN:
time: 0.001 time: 0.001 # time the SN explodes (internal units)
energy: 1.0 energy: 1.0 # energy of the explosion (internal units)
x: 0.5 x: 0.5 # x-position of explostion (internal units)
y: 0.5 y: 0.5 # y-position of explostion (internal units)
z: 0.5 z: 0.5 # z-position of explostion (internal units)
...@@ -2,6 +2,7 @@ ...@@ -2,6 +2,7 @@
# This file is part of SWIFT. # This file is part of SWIFT.
# Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk), # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@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 # 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 # it under the terms of the GNU Lesser General Public License as published
......
...@@ -1799,14 +1799,15 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) { ...@@ -1799,14 +1799,15 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
/* Cooling tasks should depend on kick and unlock sourceterms */ /* Cooling tasks should depend on kick and unlock sourceterms */
else if (t->type == task_type_cooling) { else if (t->type == task_type_cooling) {
scheduler_addunlock(sched, t->ci->kick, t); scheduler_addunlock(sched, t->ci->kick, t);
scheduler_addunlock(sched, t->ci->sourceterms, t); scheduler_addunlock(sched, t->ci->sourceterms, t);
} }
/* source terms depend on cooling if performed, else on kick. It is the last task */ /* source terms depend on cooling if performed, else on kick. It is the last
task */
else if (t->type == task_type_sourceterms) { else if (t->type == task_type_sourceterms) {
if (e->policy == engine_policy_cooling) if (e->policy == engine_policy_cooling)
scheduler_addunlock(sched, t->ci->cooling, t); scheduler_addunlock(sched, t->ci->cooling, t);
else else
scheduler_addunlock(sched, t->ci->kick, t); scheduler_addunlock(sched, t->ci->kick, t);
} }
} }
} }
......
...@@ -133,7 +133,6 @@ void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) { ...@@ -133,7 +133,6 @@ void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) {
sourceterms->supernova.z}; sourceterms->supernova.z};
const int dimen = 3; const int dimen = 3;
const double timeBase = r->e->timeBase; const double timeBase = r->e->timeBase;
const float CFL_condition = r->e->hydro_properties->CFL_condition;
TIMER_TIC; TIMER_TIC;
...@@ -145,95 +144,80 @@ void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) { ...@@ -145,95 +144,80 @@ void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) {
} }
/* is supernova still active? */ /* is supernova still active? */
if(sourceterms->supernova.status == supernova_is_not_done) if (sourceterms->supernova.status == supernova_is_not_done) {
{ /* does cell contain explosion? */
/* does cell contain explosion? */ if (count > 0) {
if (count > 0) { const int incell = is_in_cell(cell_min, cell_width, location, dimen);
const int incell = is_in_cell(cell_min, cell_width, location, dimen); if (incell == 1) {
if (incell == 1) {
/* inject SN energy into particle with highest id, if it is active */
/* inject SN energy into particle with highest id, if it is active */ int imax = 0;
int imax = 0; struct part *restrict p_sn = NULL;
struct part *restrict p_sn = NULL; struct xpart *restrict xp_sn = NULL;
struct xpart *restrict xp_sn = NULL;
for (int i = 0; i < count; i++) {
for (int i = 0; i < count; i++) {
/* Get a direct pointer on the part. */
/* Get a direct pointer on the part. */ struct part *restrict p = &parts[i];
struct part *restrict p = &parts[i]; if (p->id > imax) {
if (p->id > imax) { imax = p->id;
imax = p->id; p_sn = p;
p_sn = p; xp_sn = &xparts[i];
xp_sn = &xparts[i]; }
} }
}
/* Is this part within the time step? */
/* Is this part within the time step? */ if (p_sn->ti_begin == ti_current) {
if (p_sn->ti_begin == ti_current) {
/* Does this time step straddle the feedback injection time? */
/* Does this time step straddle the feedback injection time? */ const float t_begin = p_sn->ti_begin * timeBase;
const float t_begin = p_sn->ti_begin * timeBase; const float t_end = p_sn->ti_end * timeBase;
const float t_end = p_sn->ti_end * timeBase; if (t_begin <= sourceterms->supernova.time &&
if (t_begin <= sourceterms->supernova.time && t_end > sourceterms->supernova.time) {
t_end > sourceterms->supernova.time) {
/* store old time step */
/* add supernova feedback */ const int dti_old = p_sn->ti_end - p_sn->ti_begin;
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 */ /* add supernova feedback */
const float u_old = hydro_get_internal_energy(p_sn, 0.0); do_supernova_feedback(sourceterms, p_sn);
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;
/* label supernova as done */ message(" applied super nova, time = %d, location= %e %e %e",
sourceterms->supernova.status = supernova_is_done; ti_current, p_sn->x[0], p_sn->x[1], p_sn->x[2]);
message(" applied super nova, time = %d", ti_current); message(" applied super nova, velocity = %e %e %e", p_sn->v[0],
/* recompute time step */ p_sn->v[1], p_sn->v[2]);
const float u = hydro_get_internal_energy(p_sn, 0.0); error("end");
const float ci = sqrt(hydro_gamma_minus_one * hydro_gamma * u);
/* following needs updating to be amde consisten /* update timestep if new time step shorter than old time step */
const float ci = hydro_get_soundspeed((p_sn, 0.0); const int dti = get_part_timestep(p_sn, xp_sn, r->e);
*/ if (dti < dti_old) {
const float max_timestep = CFL_condition * p_sn->h / ci; /* new maximum timestep */ p_sn->ti_end = p_sn->ti_begin + dti;
const int old_dti = p_sn->ti_end - p_sn->ti_begin; message(" changed timestep from %d to %d", dti_old, dti);
const float old_timestep = (p_sn->ti_end - p_sn->ti_begin) * timeBase; /* current time step */
float sn_timestep = old_timestep; /* apply simple time-step limiter on all particles in same cell:
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); int i_limit = 0;
message(" old dt= %e new dt = %e", old_timestep, max_timestep); for (int i = 0; i < count; i++) {
if(max_timestep < old_timestep) struct part *restrict p = &parts[i];
{ const int dti_old = p->ti_end - p->ti_begin;
sn_timestep = max_timestep; if (dti_old > 2 * dti) {
const int new_dti = i_limit++;
get_integer_timestep(max_timestep, p_sn->ti_begin, p_sn->ti_end, r->e->timeBase_inv); const int dti_new = 2 * dti;
p_sn->ti_end = p_sn->ti_begin + new_dti; p->ti_end = p->ti_begin + dti_new;
/* apply new time-step */ message(" old step = %d new step = %d", dti_old, dti_new);
p_sn->ti_end = p_sn->ti_begin + new_dti; } else
message(" changed timestep from %d %d to %d", orig_dti, old_dti, new_dti); message(" old step = %d", dti_old);
}
/* apply simple time-step limiter on all particles in this cell: */ message(" count= %d limited timestep of %d particles ", count,
int i_limit=0; i_limit);
for (int i=0; i<count; i++) }
{ // error(" check! ");
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 (timer) TIMER_TOC(timer_dosource); if (timer) TIMER_TOC(timer_dosource);
} }
......
...@@ -36,7 +36,7 @@ ...@@ -36,7 +36,7 @@
* @param source the structure that has all the source term properties * @param source the structure that has all the source term properties
*/ */
void sourceterms_init(const struct swift_params* parameter_file, void sourceterms_init(const struct swift_params* parameter_file,
struct UnitSystem* us, struct sourceterms* source) { struct UnitSystem* us, struct sourceterms* source) {
#ifdef SN_FEEDBACK #ifdef SN_FEEDBACK
source->supernova.time = parser_get_param_double(parameter_file, "SN:time"); source->supernova.time = parser_get_param_double(parameter_file, "SN:time");
......
...@@ -31,7 +31,7 @@ struct sourceterms { ...@@ -31,7 +31,7 @@ struct sourceterms {
}; };
void sourceterms_init(const struct swift_params* parameter_file, 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); void sourceterms_print(struct sourceterms* source);
#ifdef SN_FEEDBACK #ifdef SN_FEEDBACK
#include "sourceterms/sn_feedback/sn_feedback.h" #include "sourceterms/sn_feedback/sn_feedback.h"
......
...@@ -38,7 +38,8 @@ __attribute__((always_inline)) INLINE static int is_in_cell( ...@@ -38,7 +38,8 @@ __attribute__((always_inline)) INLINE static int is_in_cell(
* @param sourceterms the structure describing the source terms properties * @param sourceterms the structure describing the source terms properties
* @param p the particle to apply feedback to * @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. * by supernova energy / particle mass.
*/ */
__attribute__((always_inline)) INLINE static void do_supernova_feedback( __attribute__((always_inline)) INLINE static void do_supernova_feedback(
...@@ -49,7 +50,9 @@ __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); hydro_set_internal_energy(p, u_new);
const float u_set = hydro_get_internal_energy(p, 0.0); const float u_set = hydro_get_internal_energy(p, 0.0);
message(" unew = %e %e s= %e", u_new, u_set, p->entropy); 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 */ #endif /* SWIFT_SN_FEEDBACK_H */
...@@ -20,21 +20,20 @@ ...@@ -20,21 +20,20 @@
* @file src/sourceterms/sn_feedback_struct.h * @file src/sourceterms/sn_feedback_struct.h
* @brief Routines related to source terms (feedback) * @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 #ifndef SWIFT_SN_FEEDBACK_STRUCT_H
#define SWIFT_SN_FEEDBACK_STRUCT_H #define SWIFT_SN_FEEDBACK_STRUCT_H
enum supernova_status { enum supernova_status { supernova_is_done, supernova_is_not_done };
supernova_is_done,
supernova_is_not_done
};
/** /**
* @file src/sourceterms/sn_feedback_struct.h * @file src/sourceterms/sn_feedback_struct.h
* @brief Routines related to source terms (feedback) * @brief Routines related to source terms (feedback)
* *
* The structure that describes the source term (supernova 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 * that records the status of the supernova
*/ */
struct supernova_struct { struct supernova_struct {
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment