Commit 45ee7266 authored by Loic Hausammann's avatar Loic Hausammann
Browse files

Start implementing star_ghost

parent 29c165dd
......@@ -96,6 +96,9 @@
/* Import the gravity loop functions. */
#include "runner_doiact_grav.h"
/* Import the star loop functions. */
#include "runner_doiact_star.h"
/**
* @brief Perform source terms
*
......@@ -132,6 +135,220 @@ void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) {
if (timer) TIMER_TOC(timer_dosource);
}
/**
* @brief Intermediate task after the density to check that the smoothing
* lengths are correct.
*
* @param r The runner thread.
* @param c The cell.
* @param timer Are we timing this ?
*/
void runner_do_star_ghost(struct runner *r, struct cell *c, int timer) {
struct spart *restrict sparts = c->sparts;
const struct engine *e = r->e;
const struct space *s = e->s;
const struct cosmology *cosmo = e->cosmology;
const float star_h_max = e->star_properties->h_max;
const float eps = e->star_properties->h_tolerance;
const float star_eta_dim =
pow_dimension(e->star_properties->eta_neighbours);
const int max_smoothing_iter = e->star_properties->max_smoothing_iterations;
int redo = 0, count = 0;
TIMER_TIC;
/* Anything to do here? */
if (!cell_is_active_star(c, e)) return;
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) runner_do_star_ghost(r, c->progeny[k], 0);
} else {
/* Init the list of active particles that have to be updated. */
int *sid = NULL;
if ((sid = (int *)malloc(sizeof(int) * c->scount)) == NULL)
error("Can't allocate memory for pid.");
for (int k = 0; k < c->scount; k++)
if (spart_is_active(&sparts[k], e)) {
sid[count] = k;
++count;
}
/* While there are particles that need to be updated... */
for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
num_reruns++) {
/* Reset the redo-count. */
redo = 0;
/* Loop over the remaining active parts in this cell. */
for (int i = 0; i < count; i++) {
/* Get a direct pointer on the part. */
struct spart *sp = &sparts[sid[i]];
#ifdef SWIFT_DEBUG_CHECKS
/* Is this part within the timestep? */
if (!spart_is_active(sp, e)) error("Ghost applied to inactive particle");
#endif
/* Get some useful values */
const float h_old = sp->h;
const float h_old_dim = pow_dimension(h_old);
const float h_old_dim_minus_one = pow_dimension_minus_one(h_old);
float h_new;
int has_no_neighbours = 0;
if (sp->wcount == 0.f) { /* No neighbours case */
/* Flag that there were no neighbours */
has_no_neighbours = 1;
/* Double h and try again */
h_new = 2.f * h_old;
} else {
/* Finish the density calculation */
star_end_density(p, cosmo);
/* Compute one step of the Newton-Raphson scheme */
const float n_sum = sp->wcount * h_old_dim;
const float n_target = star_eta_dim;
const float f = n_sum - n_target;
const float f_prime =
sp->wcount_dh * h_old_dim +
hydro_dimension * sp->wcount * h_old_dim_minus_one;
/* Avoid floating point exception from f_prime = 0 */
h_new = h_old - f / (f_prime + FLT_MIN);
#ifdef SWIFT_DEBUG_CHECKS
if ((f > 0.f && h_new > h_old) || (f < 0.f && h_new < h_old))
error(
"Smoothing length correction not going in the right direction");
#endif
/* Safety check: truncate to the range [ h_old/2 , 2h_old ]. */
h_new = min(h_new, 2.f * h_old);
h_new = max(h_new, 0.5f * h_old);
}
/* Check whether the particle has an inappropriate smoothing length */
if (fabsf(h_new - h_old) > eps * h_old) {
/* Ok, correct then */
sp->h = h_new;
/* If below the absolute maximum, try again */
if (sp->h < star_h_max) {
/* Flag for another round of fun */
sid[redo] = sid[i];
redo += 1;
/* Re-initialise everything */
star_init_spart(sp, hs);
/* Off we go ! */
continue;
} else {
/* Ok, this particle is a lost cause... */
sp->h = star_h_max;
/* Do some damage control if no neighbours at all were found */
if (has_no_neighbours) {
star_spart_has_no_neighbours(sp, cosmo);
}
}
}
/* We now have a particle whose smoothing length has converged */
/* As of here, particle force variables will be set. */
/* Compute variables required for the force loop */
star_prepare_force(sp, cosmo);
/* The particle force values are now set. Do _NOT_
try to read any particle density variables! */
/* Prepare the particle for the force loop over neighbours */
star_reset_acceleration(sp);
}
/* We now need to treat the particles whose smoothing length had not
* converged again */
/* Re-set the counter for the next loop (potentially). */
count = redo;
if (count > 0) {
/* Climb up the cell hierarchy. */
for (struct cell *finger = c; finger != NULL; finger = finger->parent) {
/* Run through this cell's density interactions. */
for (struct link *l = finger->density; l != NULL; l = l->next) {
#ifdef SWIFT_DEBUG_CHECKS
if (l->t->ti_run < r->e->ti_current)
error("Density task should have been run.");
#endif
/* Self-interaction? */
if (l->t->type == task_type_self)
runner_doself_subset_star_density(r, finger, parts, pid, count);
/* Otherwise, pair interaction? */
else if (l->t->type == task_type_pair) {
/* Left or right? */
if (l->t->ci == finger)
runner_dopair_subset_star_density(r, finger, parts, pid,
count, l->t->cj);
else
runner_dopair_subset_star_density(r, finger, parts, pid,
count, l->t->ci);
}
/* Otherwise, sub-self interaction? */
else if (l->t->type == task_type_sub_self)
runner_dosub_subset_density(r, finger, parts, pid, count, NULL,
-1, 1);
/* Otherwise, sub-pair interaction? */
else if (l->t->type == task_type_sub_pair) {
/* Left or right? */
if (l->t->ci == finger)
runner_dosub_subset_density(r, finger, parts, pid, count,
l->t->cj, -1, 1);
else
runner_dosub_subset_density(r, finger, parts, pid, count,
l->t->ci, -1, 1);
}
}
}
}
}
if (count) {
error("Smoothing length failed to converge on %i particles.", count);
}
/* Be clean */
free(sid);
}
if (timer) TIMER_TOC(timer_do_star_ghost);
}
/**
* @brief Calculate gravity acceleration from external potential
*
......@@ -2145,6 +2362,8 @@ void *runner_main(void *data) {
runner_doself_recursive_grav(r, ci, 1);
else if (t->subtype == task_subtype_external_grav)
runner_do_grav_external(r, ci, 1);
else if (t->subtype == task_subtype_star_density)
runner_doself_star_density(r, ci, 1);
else
error("Unknown/invalid task subtype (%d).", t->subtype);
break;
......@@ -2160,6 +2379,8 @@ void *runner_main(void *data) {
runner_dopair2_branch_force(r, ci, cj);
else if (t->subtype == task_subtype_grav)
runner_dopair_recursive_grav(r, ci, cj, 1);
else if (t->subtype == task_subtype_star_density)
runner_dopair_star_density(r, ci, cj, 1);
else
error("Unknown/invalid task subtype (%d).", t->subtype);
break;
......
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
* 2018 Loic Hausammann (loic.hausammann@epfl.ch)
*
* 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/>.
*
******************************************************************************/
#include "swift.h"
/**
* @brief Calculate the number density of #part around the #spart
*
* @param r runner task
* @param c cell
* @param timer 1 if the time is to be recorded.
*/
void runner_doself_star_density(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
TIMER_TIC;
/* Anything to do here? */
if (!cell_is_active_star(c, e)) return;
/* Cosmological terms */
const float a = cosmo->a;
const float H = cosmo->H;
const int scount = c->scount;
const int count = c->count;
struct spart *restrict sparts = c->sparts;
struct part *restrict parts = c->parts;
/* Loop over the sparts in ci. */
for (int sid = 0; sid < scount; sid++) {
/* Get a hold of the ith spart in ci. */
struct spart *restrict si = &sparts[sid];
const float hi = si->h;
const float hig2 = hi * hi * kernel_gamma2;
const float six[3] = {(float)(si->x[0] - c->loc[0]),
(float)(si->x[1] - c->loc[1]),
(float)(si->x[2] - c->loc[2])};
/* Loop over the parts in cj. */
for (int pjd = 0; pjd < count; pjd++) {
/* Get a pointer to the jth particle. */
struct part *restrict pj = &parts[pjd];
const float hj = pj->h;
/* Compute the pairwise distance. */
const float pjx[3] = {(float)(pj->x[0] - c->loc[0]),
(float)(pj->x[1] - c->loc[1]),
(float)(pj->x[2] - c->loc[2])};
float dx[3] = {six[0] - pjx[0], six[1] - pjx[1], six[2] - pjx[2]};
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
#ifdef SWIFT_DEBUG_CHECKS
/* Check that particles have been drifted to the current time */
if (si->ti_drift != e->ti_current)
error("Particle si not drifted to current time");
if (pj->ti_drift != e->ti_current)
error("Particle pj not drifted to current time");
#endif
if (r2 < hig2)
runner_iact_nonsym_star_density(r2, dx, hi, hj, si, pj, a, H);
} /* loop over the parts in ci. */
} /* loop over the sparts in ci. */
TIMER_TOC(timer_doself_star_density);
}
/**
* @brief Calculate the number density of #part around the #spart
*
* @param r runner task
* @param c cell
* @param timer 1 if the time is to be recorded.
*/
void runner_dosubpair_star_density(struct runner *r, struct cell *restrict ci,
struct cell *restrict cj) {
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
/* Anything to do here? */
if (!cell_is_active_star(ci, e) && !cell_is_active_star(cj, e)) return;
const int scount_i = ci->scount;
const int count_j = cj->count;
struct spart *restrict sparts_i = ci->sparts;
struct part *restrict parts_j = cj->parts;
/* Cosmological terms */
const float a = cosmo->a;
const float H = cosmo->H;
/* Get the relative distance between the pairs, wrapping. */
double shift[3] = {0.0, 0.0, 0.0};
for (int k = 0; k < 3; k++) {
if (cj->loc[k] - ci->loc[k] < -e->s->dim[k] / 2)
shift[k] = e->s->dim[k];
else if (cj->loc[k] - ci->loc[k] > e->s->dim[k] / 2)
shift[k] = -e->s->dim[k];
}
/* Loop over the sparts in ci. */
for (int sid = 0; sid < scount_i; sid++) {
/* Get a hold of the ith spart in ci. */
struct spart *restrict si = &sparts_i[sid];
const float hi = si->h;
const float hig2 = hi * hi * kernel_gamma2;
const float six[3] = {(float)(si->x[0] - (cj->loc[0] + shift[0])),
(float)(si->x[1] - (cj->loc[1] + shift[1])),
(float)(si->x[2] - (cj->loc[2] + shift[2]))};
/* Loop over the parts in cj. */
for (int pjd = 0; pjd < count_j; pjd++) {
/* Get a pointer to the jth particle. */
struct part *restrict pj = &parts_j[pjd];
const float hj = pj->h;
/* Compute the pairwise distance. */
const float pjx[3] = {(float)(pj->x[0] - cj->loc[0]),
(float)(pj->x[1] - cj->loc[1]),
(float)(pj->x[2] - cj->loc[2])};
float dx[3] = {six[0] - pjx[0], six[1] - pjx[1], six[2] - pjx[2]};
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
#ifdef SWIFT_DEBUG_CHECKS
/* Check that particles have been drifted to the current time */
if (si->ti_drift != e->ti_current)
error("Particle si not drifted to current time");
if (pj->ti_drift != e->ti_current)
error("Particle pj not drifted to current time");
#endif
if (r2 < hig2)
runner_iact_nonsym_star_density(r2, dx, hj, hi, si, pj, a, H);
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
}
void runner_dopair_star_density(struct runner *r, struct cell *restrict ci,
struct cell *restrict cj, int timer) {
TIMER_TIC;
runner_dosubpair_star_density(r, ci, cj);
runner_dosubpair_star_density(r, cj, ci);
if (timer) TIMER_TOC(timer_dopair_star_density);
}
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