Commit 2fb51c3f authored by Tom Theuns's avatar Tom Theuns
Browse files

added source terms

parent 352f3fa0
base = 'Feedback'
inf = 'Feedback_003.hdf5'
blast = [0.5,0.5,0.5] ; 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]
;; distnace 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
# 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
energy: 1.0
x: 0.5
y: 0.5
z: 0.5
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_64.hdf5
This diff is collapsed.
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@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
# 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()
...@@ -73,6 +73,8 @@ void print_help_message() { ...@@ -73,6 +73,8 @@ void print_help_message() {
"Overwrite the CPU frequency (Hz) to be used for time measurements"); "Overwrite the CPU frequency (Hz) to be used for time measurements");
printf(" %2s %8s %s\n", "-g", "", printf(" %2s %8s %s\n", "-g", "",
"Run with an external gravitational potential"); "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", "-G", "", "Run with self-gravity");
printf(" %2s %8s %s\n", "-n", "{int}", printf(" %2s %8s %s\n", "-n", "{int}",
"Execute a fixed number of time steps"); "Execute a fixed number of time steps");
...@@ -143,6 +145,7 @@ int main(int argc, char *argv[]) { ...@@ -143,6 +145,7 @@ int main(int argc, char *argv[]) {
int nsteps = -2; int nsteps = -2;
int with_cosmology = 0; int with_cosmology = 0;
int with_external_gravity = 0; int with_external_gravity = 0;
int with_sourceterms = 0;
int with_self_gravity = 0; int with_self_gravity = 0;
int with_hydro = 0; int with_hydro = 0;
int with_fp_exceptions = 0; int with_fp_exceptions = 0;
...@@ -153,7 +156,7 @@ int main(int argc, char *argv[]) { ...@@ -153,7 +156,7 @@ int main(int argc, char *argv[]) {
/* Parse the parameters */ /* Parse the parameters */
int c; int c;
while ((c = getopt(argc, argv, "acdef:gGhn:st:v:y:")) != -1) switch (c) { while ((c = getopt(argc, argv, "acdef:gFGhn:st:v:y:")) != -1) switch (c) {
case 'a': case 'a':
with_aff = 1; with_aff = 1;
break; break;
...@@ -178,6 +181,9 @@ int main(int argc, char *argv[]) { ...@@ -178,6 +181,9 @@ int main(int argc, char *argv[]) {
break; break;
case 'G': case 'G':
with_self_gravity = 1; with_self_gravity = 1;
break;
case 'F':
with_sourceterms = 1;
break; break;
case 'h': case 'h':
if (myrank == 0) print_help_message(); if (myrank == 0) print_help_message();
...@@ -335,6 +341,11 @@ int main(int argc, char *argv[]) { ...@@ -335,6 +341,11 @@ int main(int argc, char *argv[]) {
if (with_external_gravity) potential_init(params, &us, &potential); if (with_external_gravity) potential_init(params, &us, &potential);
if (with_external_gravity && myrank == 0) potential_print(&potential); if (with_external_gravity && myrank == 0) potential_print(&potential);
/* Initialise the feedback properties */
struct sourceterms sourceterms;
if (with_sourceterms) source_terms_init(params, &us, &sourceterms);
if (with_sourceterms && myrank == 0) source_terms_print(&sourceterms);
/* Read particles and space information from (GADGET) ICs */ /* Read particles and space information from (GADGET) ICs */
char ICfileName[200] = ""; char ICfileName[200] = "";
parser_get_param_string(params, "InitialConditions:file_name", ICfileName); parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
...@@ -440,6 +451,7 @@ int main(int argc, char *argv[]) { ...@@ -440,6 +451,7 @@ int main(int argc, char *argv[]) {
if (with_hydro) engine_policies |= engine_policy_hydro; if (with_hydro) engine_policies |= engine_policy_hydro;
if (with_self_gravity) engine_policies |= engine_policy_self_gravity; if (with_self_gravity) engine_policies |= engine_policy_self_gravity;
if (with_external_gravity) engine_policies |= engine_policy_external_gravity; if (with_external_gravity) engine_policies |= engine_policy_external_gravity;
if (with_sourceterms) engine_policies |= engine_policy_sourceterms;
if (with_cosmology) engine_policies |= engine_policy_cosmology; if (with_cosmology) engine_policies |= engine_policy_cosmology;
/* Initialize the engine with the space and policies. */ /* Initialize the engine with the space and policies. */
...@@ -447,7 +459,7 @@ int main(int argc, char *argv[]) { ...@@ -447,7 +459,7 @@ int main(int argc, char *argv[]) {
struct engine e; struct engine e;
engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, with_aff, engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, with_aff,
engine_policies, talking, &us, &prog_const, &hydro_properties, engine_policies, talking, &us, &prog_const, &hydro_properties,
&potential); &potential, &sourceterms);
if (myrank == 0) { if (myrank == 0) {
clocks_gettime(&toc); clocks_gettime(&toc);
message("engine_init took %.3f %s.", clocks_diff(&tic, &toc), message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
......
...@@ -42,7 +42,8 @@ endif ...@@ -42,7 +42,8 @@ endif
include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \ 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 \ 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 \ 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 physical_constants.h physical_constants_cgs.h potentials.h version.h hydro_properties.h \
sourceterms.h
# Common source files # Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
...@@ -50,7 +51,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \ ...@@ -50,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 \ 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 \ kernel_hydro.c tools.c part.c partition.c clocks.c parser.c \
physical_constants.c potentials.c hydro_properties.c \ physical_constants.c potentials.c hydro_properties.c \
runner_doiact_fft.c runner_doiact_fft.c sourceterms.c
# Include files for distribution, not installation. # 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 \ nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
...@@ -60,6 +61,7 @@ nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_h ...@@ -60,6 +61,7 @@ nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_h
gravity.h gravity_io.h \ gravity.h gravity_io.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/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 \ gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
sourceterms.h \
hydro.h hydro_io.h \ hydro.h hydro_io.h \
hydro/Minimal/hydro.h hydro/Minimal/hydro_iact.h hydro/Minimal/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 \ hydro/Minimal/hydro_debug.h hydro/Minimal/hydro_part.h \
......
...@@ -136,6 +136,9 @@ struct cell { ...@@ -136,6 +136,9 @@ struct cell {
/* Task for external gravity */ /* Task for external gravity */
struct task *grav_external; struct task *grav_external;
/* Task for source terms */
struct task *sourceterms;
/* Number of tasks that are associated with this cell. */ /* Number of tasks that are associated with this cell. */
int nr_tasks; int nr_tasks;
......
...@@ -73,6 +73,9 @@ ...@@ -73,6 +73,9 @@
#define EXTERNAL_POTENTIAL_POINTMASS #define EXTERNAL_POTENTIAL_POINTMASS
//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL //#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
/* Source terms */
#define SN_FEEDBACK
/* Are we debugging ? */ /* Are we debugging ? */
//#define SWIFT_DEBUG_CHECKS //#define SWIFT_DEBUG_CHECKS
......
...@@ -68,7 +68,7 @@ ...@@ -68,7 +68,7 @@
#include "units.h" #include "units.h"
#include "version.h" #include "version.h"
const char *engine_policy_names[13] = {"none", const char *engine_policy_names[14] = {"none",
"rand", "rand",
"steal", "steal",
"keep", "keep",
...@@ -80,7 +80,8 @@ const char *engine_policy_names[13] = {"none", ...@@ -80,7 +80,8 @@ const char *engine_policy_names[13] = {"none",
"hydro", "hydro",
"self_gravity", "self_gravity",
"external_gravity", "external_gravity",
"cosmology_integration"}; "cosmology_integration",
"sourceterms"};
/** The rank of the engine as a global variable (for messages). */ /** The rank of the engine as a global variable (for messages). */
int engine_rank; int engine_rank;
...@@ -163,6 +164,7 @@ void engine_make_gravity_hierarchical_tasks(struct engine *e, struct cell *c, ...@@ -163,6 +164,7 @@ void engine_make_gravity_hierarchical_tasks(struct engine *e, struct cell *c,
if (is_with_external_gravity) if (is_with_external_gravity)
c->grav_external = scheduler_addtask( c->grav_external = scheduler_addtask(
s, task_type_grav_external, task_subtype_none, 0, 0, c, NULL, 0); s, task_type_grav_external, task_subtype_none, 0, 0, c, NULL, 0);
} }
} }
...@@ -191,6 +193,7 @@ void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c, ...@@ -191,6 +193,7 @@ void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c,
struct scheduler *s = &e->sched; struct scheduler *s = &e->sched;
const int is_fixdt = (e->policy & engine_policy_fixdt) == engine_policy_fixdt; const int is_fixdt = (e->policy & engine_policy_fixdt) == engine_policy_fixdt;
const int is_with_sourceterms = (e->policy & engine_policy_sourceterms) == engine_policy_sourceterms;
/* Is this the super-cell? */ /* Is this the super-cell? */
if (super == NULL && (c->density != NULL || (c->count > 0 && !c->split))) { if (super == NULL && (c->density != NULL || (c->count > 0 && !c->split))) {
...@@ -225,6 +228,11 @@ void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c, ...@@ -225,6 +228,11 @@ void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c,
/* Generate the ghost task. */ /* Generate the ghost task. */
c->ghost = scheduler_addtask(s, task_type_ghost, task_subtype_none, 0, 0, c->ghost = scheduler_addtask(s, task_type_ghost, task_subtype_none, 0, 0,
c, NULL, 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);
} }
} }
...@@ -1610,6 +1618,10 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) { ...@@ -1610,6 +1618,10 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
else if (t->type == task_type_grav_external) { else if (t->type == task_type_grav_external) {
scheduler_addunlock(sched, t->ci->init, t); scheduler_addunlock(sched, t->ci->init, t);
scheduler_addunlock(sched, t, t->ci->kick); 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);
} }
} }
} }
...@@ -2602,6 +2614,11 @@ void engine_step(struct engine *e) { ...@@ -2602,6 +2614,11 @@ void engine_step(struct engine *e) {
mask |= 1 << task_type_grav_external; mask |= 1 << task_type_grav_external;
} }
/* 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 */ /* Add MPI tasks if need be */
if (e->policy & engine_policy_mpi) { if (e->policy & engine_policy_mpi) {
...@@ -2927,7 +2944,8 @@ void engine_init(struct engine *e, struct space *s, ...@@ -2927,7 +2944,8 @@ void engine_init(struct engine *e, struct space *s,
const struct UnitSystem *internal_units, const struct UnitSystem *internal_units,
const struct phys_const *physical_constants, const struct phys_const *physical_constants,
const struct hydro_props *hydro, const struct hydro_props *hydro,
const struct external_potential *potential) { const struct external_potential *potential,
const struct sourceterms *sourceterms) {
/* Clean-up everything */ /* Clean-up everything */
bzero(e, sizeof(struct engine)); bzero(e, sizeof(struct engine));
...@@ -2976,6 +2994,7 @@ void engine_init(struct engine *e, struct space *s, ...@@ -2976,6 +2994,7 @@ void engine_init(struct engine *e, struct space *s,
e->physical_constants = physical_constants; e->physical_constants = physical_constants;
e->hydro_properties = hydro; e->hydro_properties = hydro;
e->external_potential = potential; e->external_potential = potential;
e->sourceterms = sourceterms;
e->parameter_file = params; e->parameter_file = params;
engine_rank = nodeID; engine_rank = nodeID;
......
...@@ -41,6 +41,7 @@ ...@@ -41,6 +41,7 @@
#include "parser.h" #include "parser.h"
#include "partition.h" #include "partition.h"
#include "potentials.h" #include "potentials.h"
#include "sourceterms.h"
#include "runner.h" #include "runner.h"
#include "scheduler.h" #include "scheduler.h"
#include "space.h" #include "space.h"
...@@ -61,7 +62,8 @@ enum engine_policy { ...@@ -61,7 +62,8 @@ enum engine_policy {
engine_policy_hydro = (1 << 8), engine_policy_hydro = (1 << 8),
engine_policy_self_gravity = (1 << 9), engine_policy_self_gravity = (1 << 9),
engine_policy_external_gravity = (1 << 10), engine_policy_external_gravity = (1 << 10),
engine_policy_cosmology = (1 << 11) engine_policy_cosmology = (1 << 11),
engine_policy_sourceterms = (1 << 12)
}; };
extern const char *engine_policy_names[]; extern const char *engine_policy_names[];
...@@ -204,6 +206,9 @@ struct engine { ...@@ -204,6 +206,9 @@ struct engine {
/* Properties of external gravitational potential */ /* Properties of external gravitational potential */
const struct external_potential *external_potential; const struct external_potential *external_potential;
/* Properties of source terms */
const struct sourceterms *sourceterms;
/* The (parsed) parameter file */ /* The (parsed) parameter file */
const struct swift_params *parameter_file; const struct swift_params *parameter_file;
}; };
...@@ -218,7 +223,8 @@ void engine_init(struct engine *e, struct space *s, ...@@ -218,7 +223,8 @@ void engine_init(struct engine *e, struct space *s,
const struct UnitSystem *internal_units, const struct UnitSystem *internal_units,
const struct phys_const *physical_constants, const struct phys_const *physical_constants,
const struct hydro_props *hydro, const struct hydro_props *hydro,
const struct external_potential *potential); const struct external_potential *potential,
const struct sourceterms *sourceterms);
void engine_launch(struct engine *e, int nr_runners, unsigned int mask, void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
unsigned int submask); unsigned int submask);
void engine_prepare(struct engine *e); void engine_prepare(struct engine *e);
......
...@@ -40,7 +40,9 @@ ...@@ -40,7 +40,9 @@
/* ------------------------------------------------------------------------- */ /* ------------------------------------------------------------------------- */
#if defined(EOS_IDEAL_GAS) #if defined(EOS_IDEAL_GAS)
#if defined(EOS_ISOTHERMAL_GAS)
#error can't have both!
#endif
/** /**
* @brief Returns the internal energy given density and entropy * @brief Returns the internal energy given density and entropy
* *
......
...@@ -16,7 +16,8 @@ ...@@ -16,7 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. * along with this program. If not, see <http://www.gnu.org/licenses/>.
* *
******************************************************************************/ ******************************************************************************/
#ifndef SWIFT_GADGET_HYDRO_H
#define SWIFT_GADGET_HYDRO_H
#include "adiabatic_index.h" #include "adiabatic_index.h"
#include "dimension.h" #include "dimension.h"
#include "equation_of_state.h" #include "equation_of_state.h"
...@@ -361,3 +362,4 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( ...@@ -361,3 +362,4 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
/* We read u in the entropy field. We now get S from u */ /* We read u in the entropy field. We now get S from u */
p->entropy = gas_entropy_from_internal_energy(p->rho, p->entropy); p->entropy = gas_entropy_from_internal_energy(p->rho, p->entropy);
} }
#endif /* SWIFT_GADGET_HYDRO_H */
...@@ -47,6 +47,7 @@ ...@@ -47,6 +47,7 @@
#include "engine.h" #include "engine.h"
#include "error.h" #include "error.h"
#include "gravity.h" #include "gravity.h"
#include "sourceterms.h"
#include "hydro.h" #include "hydro.h"
#include "hydro_properties.h" #include "hydro_properties.h"
#include "kick.h" #include "kick.h"
...@@ -91,6 +92,70 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, ...@@ -91,6 +92,70 @@ 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_fft.h"
#include "runner_doiact_grav.h" #include "runner_doiact_grav.h"
/**
* @brief Calculate gravity acceleration from external potential
*
* @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) {
struct part *restrict parts = c->parts;
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;
const double location[3] = {sourceterms->supernova.x, sourceterms->supernova.y, sourceterms->supernova.z};
const int dimen=3;
const double timeBase = r->e->timeBase;
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;
}