Skip to content
Snippets Groups Projects
Commit df7de486 authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Use local function when recursing for drifts 

Now at least as fast as tkask version
parent e872b4b8
No related branches found
No related tags found
1 merge request!176Tasks cleanup
Hide whitespace changes
Inline Side-by-side
Showing
with 135 additions and 130 deletions
src/runner.c
+ 135
− 130
View file @ df7de486
......@@ -584,161 +584,166 @@ void runner_do_ghost(struct runner *r, struct cell *c) {
}
/**
* @brief Mapper function to drift particles and g-particles forward in time.
* @brief Drift particles and g-particles in a cell forward in time
*
* @param map_data An array of #cell%s.
* @param num_elements Chunk size.
* @param extra_data Pointer to an #engine.
* @param c The cell.
* @param e The engine.
*/
static void runner_do_drift(struct cell *c, struct engine *e) {
void runner_do_drift_mapper(void *map_data, int num_elements,
void *extra_data) {
struct engine *e = (struct engine *)extra_data;
struct cell *cells = (struct cell *)map_data;
const double timeBase = e->timeBase;
const double dt = (e->ti_current - e->ti_old) * timeBase;
const int ti_old = e->ti_old;
const int ti_current = e->ti_current;
for (int ind = 0; ind < num_elements; ind++) {
struct cell *c = &cells[ind];
if (c == NULL) continue;
/* Only drift local particles. */
if (c->nodeID != e->nodeID) continue;
struct part *const parts = c->parts;
struct xpart *const xparts = c->xparts;
struct gpart *const gparts = c->gparts;
float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
struct part *const parts = c->parts;
struct xpart *const xparts = c->xparts;
struct gpart *const gparts = c->gparts;
float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, 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};
double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, 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};
#ifdef TASK_VERBOSE
OUT;
#endif
/* No children? */
if (!c->split) {
/* Loop over all the g-particles in the cell */
const int nr_gparts = c->gcount;
for (size_t k = 0; k < nr_gparts; k++) {
/* Get a handle on the gpart. */
struct gpart *const gp = &gparts[k];
/* Drift... */
drift_gpart(gp, dt, timeBase, ti_old, ti_current);
/* 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 = fmaxf(dx2_max, dx2);
}
/* No children? */
if (!c->split) {
/* 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++) {
/* Loop over all the g-particles in the cell */
const int nr_gparts = c->gcount;
for (size_t k = 0; k < nr_gparts; k++) {
/* Get a handle on the part. */
struct part *const p = &parts[k];
struct xpart *const xp = &xparts[k];
/* Get a handle on the gpart. */
struct gpart *const gp = &gparts[k];
/* Drift... */
drift_part(p, xp, dt, timeBase, ti_old, ti_current);
/* Drift... */
drift_gpart(gp, 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 = fmaxf(dx2_max, dx2);
/* 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 = fmaxf(dx2_max, dx2);
}
/* Maximal smoothing length */
h_max = fmaxf(p->h, h_max);
/* 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++) {
/* Now collect quantities for statistics */
/* 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 = fmaxf(dx2_max, dx2);
/* Maximal smoothing length */
h_max = fmaxf(p->h, h_max);
/* 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 = p->mass;
/* 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 += 0.;
e_int += m * hydro_get_internal_energy(p, half_dt);
/* Collect entropy */
entropy += m * hydro_get_entropy(p, half_dt);
}
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 = p->mass;
/* Now, get the maximal particle motion from its square */
dx_max = sqrtf(dx2_max);
}
/* Collect mass */
mass += m;
/* Otherwise, aggregate data from children. */
else {
/* Collect momentum */
mom[0] += m * v[0];
mom[1] += m * v[1];
mom[2] += m * v[2];
/* Loop over the progeny and collect their data. */
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) {
struct cell *cp = c->progeny[k];
/* Recurse. */
runner_do_drift(cp, e);
dx_max = fmaxf(dx_max, cp->dx_max);
h_max = fmaxf(h_max, cp->h_max);
mass += cp->mass;
e_kin += cp->e_kin;
e_int += cp->e_int;
e_pot += cp->e_pot;
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];
}
}
/* 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]);
/* 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 = e_pot;
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];
}
/* Collect energies. */
e_kin += 0.5 * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
e_pot += 0.;
e_int += m * hydro_get_internal_energy(p, half_dt);
/**
* @brief Mapper function to drift particles and g-particles forward in time.
*
* @param map_data An array of #cell%s.
* @param num_elements Chunk size.
* @param extra_data Pointer to an #engine.
*/
/* Collect entropy */
entropy += m * hydro_get_entropy(p, half_dt);
}
void runner_do_drift_mapper(void *map_data, int num_elements,
void *extra_data) {
/* Now, get the maximal particle motion from its square */
dx_max = sqrtf(dx2_max);
}
struct engine *e = (struct engine *)extra_data;
struct cell *cells = (struct cell *)map_data;
/* Otherwise, aggregate data from children. */
else {
/* Loop over the progeny and collect their data. */
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) {
struct cell *cp = c->progeny[k];
/* Recurse. */
runner_do_drift_mapper(cp, 1, e);
dx_max = fmaxf(dx_max, cp->dx_max);
h_max = fmaxf(h_max, cp->h_max);
mass += cp->mass;
e_kin += cp->e_kin;
e_int += cp->e_int;
e_pot += cp->e_pot;
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];
}
}
for (int ind = 0; ind < num_elements; ind++) {
struct cell *c = &cells[ind];
/* 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 = e_pot;
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];
/* Only drift local particles. */
if (c != NULL && c->nodeID == e->nodeID)
runner_do_drift(c, e);
}
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment