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

combined with gizmo

parent c465da85
......@@ -73,8 +73,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", "-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");
......@@ -182,8 +181,8 @@ int main(int argc, char *argv[]) {
case 'G':
with_self_gravity = 1;
break;
case 'F':
with_sourceterms = 1;
case 'F':
with_sourceterms = 1;
break;
case 'h':
if (myrank == 0) print_help_message();
......
......@@ -42,8 +42,8 @@ endif
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 potentials.h version.h hydro_properties.h \
sourceterms.h
physical_constants.h physical_constants_cgs.h potentials.h version.h \
hydro_properties.h sourceterms.h threadpool.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
......@@ -51,7 +51,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 potentials.c hydro_properties.c \
runner_doiact_fft.c sourceterms.c
runner_doiact_fft.c sourceterms.c threadpool.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
......
......@@ -36,13 +36,17 @@
/* Time integration constants. */
#define const_max_u_change 0.1f
/* Thermal energy per unit mass used as a constant for the isothermal EoS */
#define const_isothermal_internal_energy 20.2615290634f
/* Dimensionality of the problem */
#define HYDRO_DIMENSION_3D
//#define HYDRO_DIMENSION_2D
//#define HYDRO_DIMENSION_3D
#define HYDRO_DIMENSION_2D
//#define HYDRO_DIMENSION_1D
/* Hydrodynamical adiabatic index. */
#define HYDRO_GAMMA_5_3
//#define HYDRO_GAMMA_7_5
//#define HYDRO_GAMMA_4_3
//#define HYDRO_GAMMA_2_1
......@@ -60,8 +64,24 @@
/* SPH variant to use */
//#define MINIMAL_SPH
#define GADGET2_SPH
//#define GADGET2_SPH
//#define DEFAULT_SPH
#define GIZMO_SPH
/* Riemann solver to use (GIZMO_SPH only) */
#define RIEMANN_SOLVER_EXACT
//#define RIEMANN_SOLVER_TRRS
//#define RIEMANN_SOLVER_HLLC
/* Type of gradients to use (GIZMO_SPH only) */
/* If no option is chosen, no gradients are used (first order scheme) */
//#define GRADIENTS_SPH
#define GRADIENTS_GIZMO
/* Types of slope limiter to use (GIZMO_SPH only) */
/* Different slope limiters can be combined */
#define SLOPE_LIMITER_PER_FACE
#define SLOPE_LIMITER_CELL_WIDE
/* Self gravity stuff. */
#define const_gravity_multipole_order 2
......@@ -70,8 +90,9 @@
#define const_gravity_eta 0.025f
/* External gravity properties */
#define EXTERNAL_POTENTIAL_POINTMASS
//#define EXTERNAL_POTENTIAL_POINTMASS
//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
#define EXTERNAL_POTENTIAL_DISK_PATCH
/* Source terms */
#define SN_FEEDBACK
......
......@@ -228,7 +228,11 @@ void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c,
/* Generate the ghost task. */
c->ghost = scheduler_addtask(s, task_type_ghost, task_subtype_none, 0, 0,
c, NULL, 0);
#ifdef EXTRA_HYDRO_LOOP
/* Generate the extra ghost task. */
c->extra_ghost = scheduler_addtask(s, task_type_extra_ghost,
task_subtype_none, 0, 0, c, NULL, 0);
#endif
/* add source terms */
if (is_with_sourceterms)
c->sourceterms = scheduler_addtask(s, task_type_sourceterms, task_subtype_none, 0, 0,
......
......@@ -41,9 +41,9 @@
#include "parser.h"
#include "partition.h"
#include "potentials.h"
#include "sourceterms.h"
#include "runner.h"
#include "scheduler.h"
#include "sourceterms.h"
#include "space.h"
#include "task.h"
#include "units.h"
......@@ -224,7 +224,7 @@ 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 sourceterms *sourceterms);
const struct sourceterms *sourceterms);
void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
unsigned int submask);
void engine_prepare(struct engine *e);
......
......@@ -618,3 +618,18 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
struct part* p) {
return 0.f;
}
/**
* @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
*
* @param p The particle
* @param u The new internal energy
*/
__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
struct part *restrict p, float u) {
p->conserved.energy = u;
}
......@@ -1226,19 +1226,6 @@ void *runner_main(void *data) {
case task_type_sourceterms:
runner_do_sourceterms(r, t->ci, 1);
break;
case task_type_part_sort:
space_do_parts_sort();
break;
case task_type_gpart_sort:
space_do_gparts_sort();
break;
case task_type_split_cell:
space_do_split(e->s, t->ci);
break;
case task_type_rewait:
scheduler_do_rewait((struct task *)t->ci, (struct task *)t->cj,
t->flags, t->rank);
break;
default:
error("Unknown task type.");
}
......
......@@ -26,11 +26,11 @@
/* Local includes. */
#include "const.h"
#include "equation_of_state.h"
#include "error.h"
#include "parser.h"
#include "part.h"
#include "physical_constants.h"
#include "equation_of_state.h"
#include "units.h"
/* This object's header. */
......@@ -44,15 +44,15 @@
* @param sourceterms: the structure that has all the source term properties
*/
void source_terms_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");
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.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");
#endif /* SN_FEEDBCK */
};
......@@ -64,13 +64,18 @@ void source_terms_init(const struct swift_params* parameter_file,
void source_terms_print(const struct sourceterms* source) {
#ifdef SN_FEEDBACK
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);
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 /* 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);
};
/* __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); */
/* }; */
......@@ -21,18 +21,17 @@
#include "./const.h"
/* So far only one model here */
struct sourceterms
{
struct sourceterms {
#ifdef SN_FEEDBACK
#include "sourceterms/sn_feedback/sn_feedback_struct.h"
#endif
};
void source_terms_init(const struct swift_params* parameter_file,
struct UnitSystem* us,
struct sourceterms* source);
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);
/* 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
......
......@@ -22,7 +22,7 @@
#include "feedback.h"
/**
* @brief Computes the feedback time-step of a given particle due to
* @brief Computes the feedback time-step of a given particle due to
* a source term
*
* This function only branches towards the potential chosen by the user.
......@@ -31,21 +31,18 @@
* @param phys_const The physical constants in internal units.
* @param g Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float
sourceterms_compute_timestep(const struct sourceterms *feedback,
const struct phys_const* const phys_const,
const struct part* const p){
__attribute__((always_inline)) INLINE static float sourceterms_compute_timestep(
const struct sourceterms* feedback,
const struct phys_const* const phys_const, const struct part* const p) {
float dt = FLT_MAX;
#ifdef SN_FEEDBACK
dt =
fmin(dt, sn_feedback_timestep(feedback, phys_const, p));
dt = fmin(dt, sn_feedback_timestep(feedback, phys_const, p));
#endif
}
__attribute__((always_inline)) INLINE static void feedback(const struct sourceterms *feedback,
const struct phys_const* const phys_const,
struct part* p)
__attribute__((always_inline)) INLINE static void feedback(
const struct sourceterms* feedback,
const struct phys_const* const phys_const, struct part* p)
#ifdef SN_FEEDBACK
sn_feedback(feedback, phys_const, p);
sn_feedback(feedback, phys_const, p);
#endif
}
#error: this file is no longer in use!
#error : this file is no longer in use!
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
......@@ -26,26 +26,26 @@
#include "../sourceterms.h"
void supernova_feedback_init(const struct swift_params* parameter_file,
struct UnitSystem* us,
struct supernova* sn){
sn.time = parser_get_param_double(parameter_file, "SN:time");
struct UnitSystem* us, struct supernova* sn) {
sn.time = parser_get_param_double(parameter_file, "SN:time");
sn.energy = parser_get_param_double(parameter_file, "SN:energy");
sn.x = parser_get_param_double(parameter_file, "SN:x");
sn.y = parser_get_param_double(parameter_file, "SN:y");
sn.z = parser_get_param_double(parameter_file, "SN:z");
sn.x = parser_get_param_double(parameter_file, "SN:x");
sn.y = parser_get_param_double(parameter_file, "SN:y");
sn.z = parser_get_param_double(parameter_file, "SN:z");
}
void supernova_feedback_print(const struct supernova* sn){
message(" Single SNe of energy= %e will explode at time= %e at location (%e,%e,%e)",
sn.energy, sn.time, sn.x, sn.y, sn.z);
void supernova_feedback_print(const struct supernova* sn) {
message(
" Single SNe of energy= %e will explode at time= %e at location "
"(%e,%e,%e)",
sn.energy, sn.time, sn.x, sn.y, sn.z);
};
__attribute__((always_inline)) INLINE static void do_supernova_feedback(const struct sourceterms* sourceterms, struct part* p)
{
};
__attribute__((always_inline)) INLINE static void do_supernova_feedback(
const struct sourceterms* sourceterms, struct part* p){};
__attribute__((always_inline)) INLINE void update_entropy(const sourceterms *sourceterms, struct part* p)
{
__attribute__((always_inline)) INLINE void update_entropy(
const sourceterms* sourceterms, struct part* p) {
/*updates the entropy of a particle due to feedback */
float u_old;
......@@ -55,10 +55,10 @@ __attribute__((always_inline)) INLINE void update_entropy(const sourceterms *sou
float rho = p->rho;
// u_old = old_entropy/(GAMMA_MINUS1) * pow(rho,GAMMA_MINUS1);
const float u_old = hydro_get_internal_energy(p,0); // dt = 0 because using current entropy
const float u_old =
hydro_get_internal_energy(p, 0); // dt = 0 because using current entropy
const float u_new = u_old + sourceterms->supernova.energy;
const float new_entropy = u_new*pow_minus_gamma_minus_one(-p>rho) * hydro_gamma_minus_one;
const float new_entropy =
u_new * pow_minus_gamma_minus_one(-p > rho) * hydro_gamma_minus_one;
p->entropy = new_entropy;
}
......@@ -24,36 +24,25 @@
#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);
// 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[], const int dimen)
{
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;
__attribute__((always_inline)) INLINE static int is_in_cell(
const double cell_min[], const double cell_width[], const double location[],
const int dimen) {
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;
};
__attribute__((always_inline)) INLINE static void do_supernova_feedback(const struct sourceterms* sourceterms, struct part* p)
{
const float density = p->rho;
const float u_old = hydro_get_internal_energy(p,0);
const float u_new = u_old + sourceterms->supernova.energy / p->mass;
const float entropy_old = gas_entropy_from_internal_energy(density, u_old);
const float entropy_new = calculate_entropy(entropy_old, density, u_old, u_new);
__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);
const float u_new = u_old + sourceterms->supernova.energy / hydro_get_mass(p);
hydro_set_internal_energy(p, u_new);
hydro_set_entropy(p,entropy_new);
message(" u_old= %e u_new= %e S_old= %e S_new= %e rho= %e", u_old, u_new, entropy_old, entropy_new, density);
message(" p->s %e", p->entropy);
};
#endif /* SWIFT_SN_FEEDBACK_H */
/* sn feedback */
struct {
double time;
double energy;
double x, y, z;
} supernova;
double time;
double energy;
double x, y, z;
} supernova;
......@@ -43,12 +43,12 @@
#include "partition.h"
#include "physical_constants.h"
#include "potentials.h"
#include "sourceterms.h"
#include "queue.h"
#include "runner.h"
#include "scheduler.h"
#include "serial_io.h"
#include "single_io.h"
#include "sourceterms.h"
#include "space.h"
#include "task.h"
#include "timers.h"
......
......@@ -48,11 +48,10 @@
/* Task type names. */
const char *taskID_names[task_type_count] = {
"none", "sort", "self", "pair", "sub_self",
"sub_pair", "init", "ghost", "drift", "kick",
"kick_fixdt", "send", "recv", "grav_gather_m", "grav_fft",
"grav_mm", "grav_up", "grav_external", "part_sort", "gpart_sort",
"split_cell", "rewait"};
"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", "grav_external", "sourceterms"};
const char *subtaskID_names[task_subtype_count] = {"none", "density", "force",
"grav", "tend"};
......@@ -111,6 +110,7 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
case task_type_sort:
case task_type_ghost:
case task_type_extra_ghost:
case task_type_sourceterms:
return task_action_part;
break;
......
......@@ -52,10 +52,6 @@ enum task_types {
task_type_grav_up,
task_type_grav_external,
task_type_sourceterms,
task_type_part_sort,
task_type_gpart_sort,
task_type_split_cell,
task_type_rewait,
task_type_count
};
......
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