Skip to content
GitLab
Commit 586c7573 authored by Tom Theuns's avatar Tom Theuns
Browse files

updated soruceterms

parent c155969a
Hide whitespace changes
Inline Side-by-side
src/engine.c
View file @ 586c7573
......@@ -1796,15 +1796,17 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
scheduler_addunlock(sched, t->ci->init, t);
scheduler_addunlock(sched, t, t->ci->kick);
}
/* source terms depend on kick (should eventually depend on cooling) */
else if (t->type == task_type_sourceterms) {
scheduler_addunlock(sched, t->ci->kick, t);
}
/* 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);
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);
}
}
}
......@@ -3133,7 +3135,7 @@ void engine_init(struct engine *e, struct space *s,
const struct hydro_props *hydro,
const struct external_potential *potential,
const struct cooling_data *cooling,
const struct sourceterms *sourceterms) {
struct sourceterms *sourceterms) {
/* Clean-up everything */
bzero(e, sizeof(struct engine));
......
src/engine.h
View file @ 586c7573
......@@ -210,7 +210,7 @@ struct engine {
const struct cooling_data *cooling_data;
/* Properties of source terms */
const struct sourceterms *sourceterms;
struct sourceterms *sourceterms;
/* The (parsed) parameter file */
const struct swift_params *parameter_file;
......@@ -228,7 +228,7 @@ void engine_init(struct engine *e, struct space *s,
const struct hydro_props *hydro,
const struct external_potential *potential,
const struct cooling_data *cooling,
const struct sourceterms *sourceterms);
struct sourceterms *sourceterms);
void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
unsigned int submask);
void engine_prepare(struct engine *e);
......
src/runner.c
View file @ 586c7573
......@@ -113,7 +113,7 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
#include "runner_doiact_grav.h"
/**
* @brief Calculate gravity acceleration from external potential
* @brief Perform source terms
*
* @param r runner task
* @param c cell
......@@ -122,16 +122,18 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) {
struct part *restrict parts = c->parts;
struct xpart *restrict xparts = c->xparts;
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]};
const int ti_current = r->e->ti_current;
const struct sourceterms *sourceterms = r->e->sourceterms;
struct sourceterms *sourceterms = r->e->sourceterms;
const double location[3] = {sourceterms->supernova.x,
sourceterms->supernova.y,
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;
......@@ -142,37 +144,96 @@ void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) {
return;
}
/* 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 sn_part = 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;
sn_part = p;
}
}
/* Is this part within the time step? */
if (sn_part->ti_begin == ti_current) {
/* Does this time step straddle the feedback injection time? */
const float t_begin = sn_part->ti_begin * timeBase;
const float t_end = sn_part->ti_end * timeBase;
if (t_begin <= sourceterms->supernova.time &&
t_end > sourceterms->supernova.time) {
do_supernova_feedback(sourceterms, sn_part);
}
}
}
}
/* 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 (timer) TIMER_TOC(timer_dosource);
}
......
src/sourceterms.c
View file @ 586c7573
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
* 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
......@@ -21,16 +20,9 @@
/* Config parameters. */
#include "../config.h"
/* Some standard headers. */
#include <math.h>
/* Local includes. */
#include "const.h"
#include "equation_of_state.h"
#include "error.h"
#include "parser.h"
#include "part.h"
#include "physical_constants.h"
#include "units.h"
/* This object's header. */
......@@ -53,6 +45,7 @@ void source_terms_init(const struct swift_params* parameter_file,
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;
#endif /* SN_FEEDBCK */
};
......@@ -71,11 +64,3 @@ void source_terms_print(const struct sourceterms* source) {
source->supernova.y, source->supernova.z);
#endif /* SN_FEEDBACK */
};
/* __attribute__((always_inline)) INLINE float calculate_entropy( */
/* const float entropy_old, const float density, const float u_old, */
/* const float u_new) { */
/* return entropy_old + */
/* hydro_gamma_minus_one * (u_new - u_old) * */
/* pow_minus_gamma_minus_one(density); */
/* }; */
src/sourceterms.h
View file @ 586c7573
......@@ -18,20 +18,21 @@
******************************************************************************/
#ifndef SWIFT_SOURCETERMS_H
#define SWIFT_SOURCETERMS_H
#include "./const.h"
#ifdef SN_FEEDBACK
#include "sourceterms/sn_feedback/sn_feedback_struct.h"
#endif
/* So far only one model here */
struct sourceterms {
#ifdef SN_FEEDBACK
#include "sourceterms/sn_feedback/sn_feedback_struct.h"
struct supernova_struct supernova;
#endif
};
void source_terms_init(const struct swift_params* parameter_file,
struct UnitSystem* us, struct sourceterms* source);
void source_terms_print(const struct sourceterms* source);
/* float calculate_entropy(const float entropy_old, const float density, */
/* const float u_old, const float u_new); */
#ifdef SN_FEEDBACK
#include "sourceterms/sn_feedback/sn_feedback.h"
#endif
......
src/sourceterms/sn_feedback/sn_feedback.h
View file @ 586c7573
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
* Matthieu Schaller (matthieu.schaller@durham.ac.uk)
* Richard Bower (r.g.bower@durham.ac.uk)
* Stefan Arridge (stefan.arridge@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
......@@ -24,9 +21,6 @@
#include <float.h>
#include "equation_of_state.h"
#include "hydro.h"
// void do_supernova_feedback(const struct source_terms* source_terms, const
// struct phys_const* const constants, const struct part* p);
/* determine whether location is in cell */
__attribute__((always_inline)) INLINE static int is_in_cell(
const double cell_min[], const double cell_width[], const double location[],
......@@ -38,10 +32,24 @@ __attribute__((always_inline)) INLINE static int is_in_cell(
return 1;
};
/**
* @file src/sourceterms/sn_feedback.h
* @brief Routines related to source terms (supernova feedback)
* @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
* by supernova energy / particle mass.
*/
__attribute__((always_inline)) INLINE static void do_supernova_feedback(
const struct sourceterms* sourceterms, struct part* p) {
const float u_old = hydro_get_internal_energy(p, 0);
message(" u_old= %e entropy= %e", u_old, p->entropy);
const float u_new = u_old + sourceterms->supernova.energy / hydro_get_mass(p);
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);
};
#endif /* SWIFT_SN_FEEDBACK_H */
src/sourceterms/sn_feedback/sn_feedback_struct.h
View file @ 586c7573
/* sn feedback */
struct {
/*******************************************************************************
* 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;
} supernova;
enum supernova_status status;
};
#endif SWIFT_SN_FEEDBACK_STRUCT_H
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