Commit 2f9ebf13 authored by Stefan Arridge's avatar Stefan Arridge
Browse files

Added external potential energy calculation to statistics.c

parents 83b901b1 4128aff0
......@@ -120,6 +120,10 @@ int main(int argc, char *argv[]) {
error("MPI_Comm_size failed with error %i.", res);
if ((res = MPI_Comm_rank(MPI_COMM_WORLD, &myrank)) != MPI_SUCCESS)
error("Call to MPI_Comm_rank failed with error %i.", res);
/* Make sure messages are stamped with the correct rank. */
engine_rank = myrank;
if ((res = MPI_Comm_set_errhandler(MPI_COMM_WORLD, MPI_ERRORS_RETURN)) !=
MPI_SUCCESS)
error("Call to MPI_Comm_set_errhandler failed with error %i.", res);
......@@ -131,6 +135,7 @@ int main(int argc, char *argv[]) {
message("WARNING: you should use the non-MPI version of this program.");
}
fflush(stdout);
#endif
/* Let's pin the main thread */
......@@ -572,8 +577,8 @@ int main(int argc, char *argv[]) {
fprintf(file_thread, " %03i 0 0 0 0 %lli %lli 0 0 0 0 %lli\n", myrank,
e.tic_step, e.toc_step, cpufreq);
int count = 0;
for (int l = 0; l < e.sched.nr_tasks; l++)
if (!e.sched.tasks[l].skip && !e.sched.tasks[l].implicit) {
for (int l = 0; l < e.sched.nr_tasks; l++) {
if (!e.sched.tasks[l].implicit && e.sched.tasks[l].toc != 0) {
fprintf(
file_thread, " %03i %i %i %i %i %lli %lli %i %i %i %i %i\n",
myrank, e.sched.tasks[l].rid, e.sched.tasks[l].type,
......@@ -588,11 +593,10 @@ int main(int argc, char *argv[]) {
(e.sched.tasks[l].cj != NULL) ? e.sched.tasks[l].cj->gcount
: 0,
e.sched.tasks[l].flags);
fflush(stdout);
count++;
}
message("rank %d counted %d tasks", myrank, count);
fflush(stdout);
count++;
}
fclose(file_thread);
}
......@@ -608,8 +612,8 @@ int main(int argc, char *argv[]) {
/* Add some information to help with the plots */
fprintf(file_thread, " %i %i %i %i %lli %lli %i %i %i %lli\n", -2, -1, -1,
1, e.tic_step, e.toc_step, 0, 0, 0, cpufreq);
for (int l = 0; l < e.sched.nr_tasks; l++)
if (!e.sched.tasks[l].skip && !e.sched.tasks[l].implicit)
for (int l = 0; l < e.sched.nr_tasks; l++) {
if (!e.sched.tasks[l].implicit && e.sched.tasks[l].toc != 0) {
fprintf(
file_thread, " %i %i %i %i %lli %lli %i %i %i %i\n",
e.sched.tasks[l].rid, e.sched.tasks[l].type,
......@@ -619,6 +623,8 @@ int main(int argc, char *argv[]) {
(e.sched.tasks[l].cj == NULL) ? 0 : e.sched.tasks[l].cj->count,
(e.sched.tasks[l].ci == NULL) ? 0 : e.sched.tasks[l].ci->gcount,
(e.sched.tasks[l].cj == NULL) ? 0 : e.sched.tasks[l].cj->gcount);
}
}
fclose(file_thread);
#endif
}
......
......@@ -44,7 +44,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.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 potential.h version.h \
hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h \
sourceterms_struct.h
sourceterms_struct.h statistics.h
# Common source files
......@@ -53,7 +53,8 @@ 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 potential.c hydro_properties.c \
runner_doiact_fft.c threadpool.c cooling.c sourceterms.c
runner_doiact_fft.c threadpool.c cooling.c sourceterms.c \
statistics.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
......
......@@ -52,6 +52,7 @@
#include "gravity.h"
#include "hydro.h"
#include "hydro_properties.h"
#include "scheduler.h"
#include "space.h"
#include "timers.h"
......@@ -763,6 +764,9 @@ void cell_clean_links(struct cell *c, void *data) {
c->force = NULL;
c->nr_force = 0;
c->grav = NULL;
c->nr_grav = 0;
}
/**
......@@ -893,6 +897,120 @@ int cell_is_drift_needed(struct cell *c, int ti_current) {
return 0;
}
/**
* @brief Un-skips all the tasks associated with a given cell and checks
* if the space needs to be rebuilt.
*
* @param c the #cell.
* @param s the #scheduler.
*
* @return 1 If the space needs rebuilding. 0 otherwise.
*/
int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
/* Un-skip the density tasks involved with this cell. */
for (struct link *l = c->density; l != NULL; l = l->next) {
struct task *t = l->t;
const struct cell *ci = t->ci;
const struct cell *cj = t->cj;
scheduler_activate(s, t);
/* Set the correct sorting flags */
if (t->type == task_type_pair) {
if (!(ci->sorted & (1 << t->flags))) {
atomic_or(&ci->sorts->flags, (1 << t->flags));
scheduler_activate(s, ci->sorts);
}
if (!(cj->sorted & (1 << t->flags))) {
atomic_or(&cj->sorts->flags, (1 << t->flags));
scheduler_activate(s, cj->sorts);
}
}
/* Check whether there was too much particle motion */
if (t->type == task_type_pair || t->type == task_type_sub_pair) {
if (t->tight &&
(max(ci->h_max, cj->h_max) + ci->dx_max + cj->dx_max > cj->dmin ||
ci->dx_max > space_maxreldx * ci->h_max ||
cj->dx_max > space_maxreldx * cj->h_max))
return 1;
#ifdef WITH_MPI
/* Activate the send/recv flags. */
if (ci->nodeID != engine_rank) {
/* Activate the tasks to recv foreign cell ci's data. */
scheduler_activate(s, ci->recv_xv);
scheduler_activate(s, ci->recv_rho);
scheduler_activate(s, ci->recv_ti);
/* Look for the local cell cj's send tasks. */
struct link *l = NULL;
for (l = cj->send_xv; l != NULL && l->t->cj->nodeID != ci->nodeID;
l = l->next)
;
if (l == NULL) error("Missing link to send_xv task.");
scheduler_activate(s, l->t);
for (l = cj->send_rho; l != NULL && l->t->cj->nodeID != ci->nodeID;
l = l->next)
;
if (l == NULL) error("Missing link to send_rho task.");
scheduler_activate(s, l->t);
for (l = cj->send_ti; l != NULL && l->t->cj->nodeID != ci->nodeID;
l = l->next)
;
if (l == NULL) error("Missing link to send_ti task.");
scheduler_activate(s, l->t);
} else if (cj->nodeID != engine_rank) {
/* Activate the tasks to recv foreign cell cj's data. */
scheduler_activate(s, cj->recv_xv);
scheduler_activate(s, cj->recv_rho);
scheduler_activate(s, cj->recv_ti);
/* Look for the local cell ci's send tasks. */
struct link *l = NULL;
for (l = ci->send_xv; l != NULL && l->t->cj->nodeID != cj->nodeID;
l = l->next)
;
if (l == NULL) error("Missing link to send_xv task.");
scheduler_activate(s, l->t);
for (l = ci->send_rho; l != NULL && l->t->cj->nodeID != cj->nodeID;
l = l->next)
;
if (l == NULL) error("Missing link to send_rho task.");
scheduler_activate(s, l->t);
for (l = ci->send_ti; l != NULL && l->t->cj->nodeID != cj->nodeID;
l = l->next)
;
if (l == NULL) error("Missing link to send_ti task.");
scheduler_activate(s, l->t);
}
#endif
}
}
/* Unskip all the other task types. */
for (struct link *l = c->gradient; l != NULL; l = l->next)
scheduler_activate(s, l->t);
for (struct link *l = c->force; l != NULL; l = l->next)
scheduler_activate(s, l->t);
for (struct link *l = c->grav; l != NULL; l = l->next)
scheduler_activate(s, l->t);
if (c->extra_ghost != NULL) scheduler_activate(s, c->extra_ghost);
if (c->ghost != NULL) scheduler_activate(s, c->ghost);
if (c->init != NULL) scheduler_activate(s, c->init);
if (c->kick != NULL) scheduler_activate(s, c->kick);
if (c->cooling != NULL) scheduler_activate(s, c->cooling);
if (c->sourceterms != NULL) scheduler_activate(s, c->sourceterms);
return 0;
}
/**
* @brief Set the super-cell pointers for all cells in a hierarchy.
*
......
......@@ -39,6 +39,7 @@
/* Avoid cyclic inclusions */
struct space;
struct scheduler;
/* Max tag size set to 2^29 to take into account some MPI implementations
* that use 2^31 as the upper bound on MPI tags and the fact that
......@@ -158,9 +159,6 @@ struct cell {
/* Tasks for gravity tree. */
struct task *grav_up, *grav_down;
/* Task for external gravity */
struct task *grav_external;
/* Task for cooling */
struct task *cooling;
......@@ -179,12 +177,6 @@ struct cell {
/* ID of the previous owner, e.g. runner. */
int owner;
/* Momentum of particles in cell. */
double mom[3], ang_mom[3];
/* Mass, potential, internal and kinetic energy of particles in this cell. */
double mass, e_pot_self, e_pot_ext, e_pot, e_int, e_kin, e_rad, entropy;
/* Number of particles updated in this cell. */
int updated, g_updated;
......@@ -236,6 +228,7 @@ int cell_are_neighbours(const struct cell *restrict ci,
void cell_check_multipole(struct cell *c, void *data);
void cell_clean(struct cell *c);
int cell_is_drift_needed(struct cell *c, int ti_current);
int cell_unskip_tasks(struct cell *c, struct scheduler *s);
void cell_set_super(struct cell *c, struct cell *super);
#endif /* SWIFT_CELL_H */
This diff is collapsed.
......@@ -163,7 +163,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
__attribute__((always_inline)) INLINE static float external_gravity_get_potential_energy(
const struct external_potential* potential,
const struct phys_const* const phys_const, struct part* p) {
const struct phys_const* const phys_const, const struct part* p) {
return 0.;
}
......
......@@ -122,16 +122,16 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
*
* @param potential The #external_potential used in the run.
* @param phys_const Physical constants in internal units.
* @param p Pointer to the particle data.
* @param g Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float external_gravity_get_potential_energy(
const struct external_potential* potential,
const struct phys_const* const phys_const, struct part* p) {
const struct phys_const* const phys_const, const struct gpart* g) {
const float dx = p->x[0] - potential->x;
const float dy = p->x[1] - potential->y;
const float dz = p->x[2] - potential->z;
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
return potential->vrot * potential->vrot * 0.5 * log(dx*dx + dy*dy * dz*dz);
}
......
......@@ -68,6 +68,21 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
double time, const struct external_potential* restrict potential,
const struct phys_const* restrict phys_const, struct gpart* restrict g) {}
/**
* @brief Computes the gravitational potential energy due to nothing.
*
* @param potential The #external_potential used in the run.
* @param phys_const Physical constants in internal units.
* @param g Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float external_gravity_get_potential_energy(
const struct external_potential* potential,
const struct phys_const* const phys_const, const struct part* g) {
return 0.;
}
/**
* @brief Initialises the external potential properties in the internal system
* of units.
......
......@@ -118,20 +118,20 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
/**
* @brief Computes the gravitational potential energy of a particle in an isothermal potential.
* @brief Computes the gravitational potential energy of a particle in a point mass potential.
*
* @param potential The #external_potential used in the run.
* @param phys_const Physical constants in internal units.
* @param p Pointer to the particle data.
* @param g Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float external_gravity_get_potential_energy(
const struct external_potential* potential,
const struct phys_const* const phys_const, struct part* p) {
const struct phys_const* const phys_const, const struct part* g) {
const float dx = p->x[0] - potential->x;
const float dy = p->x[1] - potential->y;
const float dz = p->x[2] - potential->z;
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
const float rinv = 1./sqrtf(dx * dx + dy * dy + dz * dz);
return -phys_const->const_newton_G * potential->mass * r_inv;
}
......
......@@ -125,16 +125,16 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
*
* @param potential The #external_potential used in the run.
* @param phys_const Physical constants in internal units.
* @param p Pointer to the particle data.
* @param g Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float external_gravity_get_potential_energy(
const struct external_potential* potential,
const struct phys_const* const phys_const, struct part* p) {
const struct phys_const* const phys_const, const struct part* g) {
const float dx = p->x[0] - potential->x;
const float dy = p->x[1] - potential->y;
const float dz = p->x[2] - potential->z;
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
return potential->vrot * potential->vrot * 0.5 * log(dx*dx + dy*dy * dz*dz + potential->epsilon2);
}
......
......@@ -748,114 +748,99 @@ void runner_do_ghost(struct runner *r, struct cell *c) {
}
/**
* @brief Drift particles and g-particles in a cell forward in time
* @brief Drift particles and g-particles in a cell forward in time,
* unskipping any tasks associated with active cells.
*
* @param c The cell.
* @param e The engine.
* @param drift whether to actually drift the particles, will not be
* necessary for non-local cells.
*/
static void runner_do_drift(struct cell *c, struct engine *e) {
static void runner_do_drift(struct cell *c, struct engine *e, int drift) {
const int ti_current = e->ti_current;
/* Unskip any active tasks. */
if (c->ti_end_min == e->ti_current) {
const int forcerebuild = cell_unskip_tasks(c, &e->sched);
if (forcerebuild) atomic_inc(&e->forcerebuild);
}
/* Do we really need to drift? */
if (drift) {
if (!e->drift_all && !cell_is_drift_needed(c, ti_current)) return;
} else {
/* Not drifting, but may still need to recurse for task skipping. */
if (c->split) {
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL) {
struct cell *cp = c->progeny[k];
runner_do_drift(cp, e, 0);
}
}
}
return;
}
const double timeBase = e->timeBase;
const int ti_old = c->ti_old;
const int ti_current = e->ti_current;
struct part *const parts = c->parts;
struct xpart *const xparts = c->xparts;
struct gpart *const gparts = c->gparts;
const struct phys_const *constants = e->physical_constants;
const struct external_potential *potential = e->external_potential;
/* Do we need to drift ? */
if (!e->drift_all && !cell_is_drift_needed(c, ti_current)) return;
/* Check that we are actually going to move forward. */
if (ti_current == ti_old) return;
/* Drift from the last time the cell was drifted to the current time */
const double dt = (ti_current - ti_old) * timeBase;
float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
double e_kin = 0.0, e_int = 0.0, e_pot_self = 0.0, e_pot_ext = 0.0, e_rad = 0.0;
double entropy = 0.0, mass = 0.0;
double mom[3] = {0.0, 0.0, 0.0};
double ang_mom[3] = {0.0, 0.0, 0.0};
/* No children? */
if (!c->split) {
/* Loop over all the g-particles in the cell */
const size_t nr_gparts = c->gcount;
for (size_t k = 0; k < nr_gparts; k++) {
/* Check that we are actually going to move forward. */
if (ti_current >= ti_old) {
/* Get a handle on the gpart. */
struct gpart *const gp = &gparts[k];
/* Loop over all the g-particles in the cell */
const size_t nr_gparts = c->gcount;
for (size_t k = 0; k < nr_gparts; k++) {
/* Drift... */
drift_gpart(gp, dt, timeBase, ti_old, ti_current);
/* Get a handle on the gpart. */
struct gpart *const gp = &gparts[k];
/* Compute (square of) motion since last cell construction */
const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
gp->x_diff[1] * gp->x_diff[1] +
gp->x_diff[2] * gp->x_diff[2];
dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
}
/* Drift... */
drift_gpart(gp, dt, timeBase, ti_old, ti_current);
/* Loop over all the particles in the cell (more work for these !) */
const size_t nr_parts = c->count;
for (size_t k = 0; k < nr_parts; k++) {
/* Compute (square of) motion since last cell construction */
const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
gp->x_diff[1] * gp->x_diff[1] +
gp->x_diff[2] * gp->x_diff[2];
dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
}
/* Get a handle on the part. */
struct part *const p = &parts[k];
struct xpart *const xp = &xparts[k];
/* Drift... */
drift_part(p, xp, dt, timeBase, ti_old, ti_current);
/* Compute (square of) motion since last cell construction */
const float dx2 = xp->x_diff[0] * xp->x_diff[0] +
xp->x_diff[1] * xp->x_diff[1] +
xp->x_diff[2] * xp->x_diff[2];
dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
/* Maximal smoothing length */
h_max = (h_max > p->h) ? h_max : p->h;
/* Now collect quantities for statistics */
const float half_dt =
(ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
const double x[3] = {p->x[0], p->x[1], p->x[2]};
const float v[3] = {xp->v_full[0] + p->a_hydro[0] * half_dt,
xp->v_full[1] + p->a_hydro[1] * half_dt,
xp->v_full[2] + p->a_hydro[2] * half_dt};
const float m = hydro_get_mass(p);
/* Collect mass */
mass += m;
/* Collect momentum */
mom[0] += m * v[0];
mom[1] += m * v[1];
mom[2] += m * v[2];
/* Collect angular momentum */
ang_mom[0] += m * (x[1] * v[2] - x[2] * v[1]);
ang_mom[1] += m * (x[2] * v[0] - x[0] * v[2]);
ang_mom[2] += m * (x[0] * v[1] - x[1] * v[0]);
/* Collect energies. */
e_kin += 0.5 * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
e_pot_self += 0.;
e_pot_ext += m * external_gravity_get_potential_energy(potential,constants,p);
e_int += m * hydro_get_internal_energy(p, half_dt);
e_rad += cooling_get_radiated_energy(xp);
/* Collect entropy */
entropy += m * hydro_get_entropy(p, half_dt);
}
/* Loop over all the particles in the cell */
const size_t nr_parts = c->count;
for (size_t k = 0; k < nr_parts; k++) {
/* Get a handle on the part. */
struct part *const p = &parts[k];
struct xpart *const xp = &xparts[k];
/* Drift... */
drift_part(p, xp, dt, timeBase, ti_old, ti_current);
/* Compute (square of) motion since last cell construction */
const float dx2 = xp->x_diff[0] * xp->x_diff[0] +
xp->x_diff[1] * xp->x_diff[1] +
xp->x_diff[2] * xp->x_diff[2];
dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
/* Maximal smoothing length */
h_max = (h_max > p->h) ? h_max : p->h;
}
/* Now, get the maximal particle motion from its square */
dx_max = sqrtf(dx2_max);
/* Now, get the maximal particle motion from its square */
dx_max = sqrtf(dx2_max);
} /* Check that we are actually going to move forward. */
}
/* Otherwise, aggregate data from children. */
......@@ -867,42 +852,15 @@ static void runner_do_drift(struct cell *c, struct engine *e) {
struct cell *cp = c->progeny[k];
/* Recurse. */
runner_do_drift(cp, e);
runner_do_drift(cp, e, drift);
dx_max = max(dx_max, cp->dx_max);
h_max = max(h_max, cp->h_max);
mass += cp->mass;
e_kin += cp->e_kin;
e_int += cp->e_int;
e_pot_self += cp->e_pot_self;
e_pot_ext += cp->e_pot_ext;
e_rad += cp->e_rad;
entropy += cp->entropy;
mom[0] += cp->mom[0];
mom[1] += cp->mom[1];
mom[2] += cp->mom[2];
ang_mom[0] += cp->ang_mom[0];
ang_mom[1] += cp->ang_mom[1];
ang_mom[2] += cp->ang_mom[2];
}
}
/* Store the values */
c->h_max = h_max;
c->dx_max = dx_max;
c->mass = mass;
c->e_kin = e_kin;
c->e_int = e_int;
c->e_pot_self = e_pot_self;
c->e_pot_ext = e_pot_ext;
c->e_rad = e_rad;
c->entropy = entropy;
c->mom[0] = mom[0];
c->mom[1] = mom[1];
c->mom[2] = mom[2];
c->ang_mom[0] = ang_mom[0];
c->ang_mom[1] = ang_mom[1];
c->ang_mom[2] = ang_mom[2];
/* Update the time of the last drift */
c->ti_old = ti_current;
......@@ -924,9 +882,11 @@ void runner_do_drift_mapper(void *map_data, int num_elements,
for (int ind = 0; ind < num_elements; ind++) {
struct cell *c = &cells[ind];
/* Only drift local particles. */
if (c != NULL && c->nodeID == e->nodeID) runner_do_drift(c, e);