Commit 0ac7431e authored by Matthieu Schaller's avatar Matthieu Schaller Committed by Peter W. Draper
Browse files

Particle splitting "a la EAGLE"

parent 8dafa3fb
......@@ -176,6 +176,18 @@ __attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
return 0.f;
}
/**
* @brief Split the coolong content of a particle into n pieces
*
* Nothing to do here.
*
* @param p The #part.
* @param xp The #xpart.
* @param n The number of pieces to split into.
*/
static INLINE void cooling_split_part(struct part* p, struct xpart* xp,
double n) {}
/**
* @brief Initialises the cooling properties.
*
......
......@@ -1629,6 +1629,17 @@ void engine_prepare(struct engine *e) {
engine_fof(e, /*dump_results=*/0, /*seed_black_holes=*/1);
}
/* Perform particle splitting. Only if there is a rebuild coming and no
repartitioning. */
if (e->forcerebuild && !e->forcerepart && e->step > 1) {
/* Let's start by drifting everybody to the current time */
if (!drifted_all) engine_drift_all(e, /*drift_mpole=*/0);
drifted_all = 1;
engine_split_gas_particles(e);
}
/* Do we need repartitioning ? */
if (e->forcerepart) {
......
......@@ -561,6 +561,9 @@ void engine_make_fof_tasks(struct engine *e);
/* Function prototypes, engine_marktasks.c. */
int engine_marktasks(struct engine *e);
/* Function prototypes, engine_split_particles.c. */
void engine_split_gas_particles(struct engine *e);
#ifdef HAVE_SETAFFINITY
cpu_set_t *engine_entry_affinity(void);
#endif
......
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl)
*
* 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/>.
*
******************************************************************************/
/* Config parameters. */
#include "../config.h"
/* This object's header. */
#include "engine.h"
/* Local headers */
#include "active.h"
#include "black_holes_struct.h"
#include "chemistry.h"
#include "cooling.h"
#include "hydro.h"
#include "random.h"
#include "star_formation.h"
#include "tracers.h"
const int particle_split_factor = 2;
const double displacement_factor = 0.2;
/**
* @brief Identify all the gas particles above a given mass threshold
* and split them into 2.
*
* This may reallocate the space arrays as new particles are created.
* This is an expensive function as it has to loop multiple times over
* the local array of particles. In case of reallocations, it may
* also have to loop over the gravity and other arrays.
*
* @param e The #engine.
*/
void engine_split_gas_particles(struct engine *e) {
/* Abort if we are not doing any splitting */
if (!e->hydro_properties->particle_splitting) return;
/* Time this */
const ticks tic = getticks();
/* Get useful constants */
struct space *s = e->s;
const int with_gravity = (e->policy & engine_policy_self_gravity) ||
(e->policy & engine_policy_external_gravity);
const double mass_threshold =
e->hydro_properties->particle_splitting_mass_threshold;
/* Quick check to avoid problems */
if (particle_split_factor != 2) {
error(
"Invalid splitting factor. Can currently only split particles into 2!");
}
/* Start by counting how many particles are above the threshold for splitting
*/
size_t counter = 0;
const size_t nr_parts_old = s->nr_parts;
const struct part *parts_old = s->parts;
for (size_t i = 0; i < nr_parts_old; ++i) {
const struct part *p = &parts_old[i];
/* Ignore inhibited particles */
if (part_is_inhibited(p, e)) continue;
/* Is the mass of this particle larger than the threshold? */
const float gas_mass = hydro_get_mass(p);
if (gas_mass > mass_threshold) ++counter;
}
/* Early abort? */
if (counter == 0) {
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
clocks_getunit());
return;
}
/* Be verbose about this. This is an important event */
message("Splitting %zd particles above the mass threshold", counter);
/* Number of particles to create */
const size_t count_new_gas = counter * particle_split_factor;
/* Do we need to reallocate the gas array for the new particles? */
if (s->nr_parts + count_new_gas > s->size_parts) {
const size_t nr_parts_new = s->nr_parts + count_new_gas;
s->size_parts = engine_parts_size_grow * nr_parts_new;
if (e->verbose) message("Reallocating the part array!");
struct part *parts_new = NULL;
if (swift_memalign("parts", (void **)&parts_new, part_align,
sizeof(struct part) * s->size_parts) != 0)
error("Failed to allocate new part data.");
memcpy(parts_new, s->parts, sizeof(struct part) * s->nr_parts);
swift_free("parts", s->parts);
struct xpart *xparts_new = NULL;
if (swift_memalign("xparts", (void **)&xparts_new, xpart_align,
sizeof(struct xpart) * s->size_parts) != 0)
error("Failed to allocate new part data.");
memcpy(xparts_new, s->xparts, sizeof(struct xpart) * s->nr_parts);
swift_free("xparts", s->xparts);
s->xparts = xparts_new;
s->parts = parts_new;
}
/* Do we need to reallocate the gpart array for the new particles? */
if (with_gravity && s->nr_gparts + count_new_gas > s->size_gparts) {
const size_t nr_gparts_new = s->nr_gparts + count_new_gas;
s->size_gparts = engine_parts_size_grow * nr_gparts_new;
if (e->verbose) message("Reallocating the gpart array!");
struct gpart *gparts_new = NULL;
if (swift_memalign("gparts", (void **)&gparts_new, gpart_align,
sizeof(struct gpart) * s->size_gparts) != 0)
error("Failed to allocate new gpart data.");
/* Offset of the new array */
const ptrdiff_t offset = gparts_new - s->gparts;
/* Copy the particles */
memcpy(gparts_new, s->gparts, sizeof(struct gpart) * s->nr_gparts);
swift_free("gparts", s->gparts);
/* We now need to correct all the pointers */
for (size_t i = 0; i < s->nr_parts; ++i) {
if (s->parts[i].time_bin <= num_time_bins) s->parts[i].gpart += offset;
}
/* We now need to correct all the pointers */
for (size_t i = 0; i < s->nr_sparts; ++i) {
if (s->sparts[i].time_bin <= num_time_bins) s->sparts[i].gpart += offset;
}
/* We now need to correct all the pointers */
for (size_t i = 0; i < s->nr_bparts; ++i) {
if (s->bparts[i].time_bin <= num_time_bins) s->bparts[i].gpart += offset;
}
s->gparts = gparts_new;
}
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that whatever reallocation happened we are still having correct
* links */
part_verify_links(s->parts, s->gparts, s->sparts, s->bparts, s->nr_parts,
s->nr_gparts, s->nr_sparts, s->nr_bparts, e->verbose);
#endif
size_t k_parts = s->nr_parts;
size_t k_gparts = s->nr_gparts;
/* Loop over the particles again to split them */
struct part *parts = s->parts;
struct xpart *xparts = s->xparts;
struct gpart *gparts = s->gparts;
for (size_t i = 0; i < nr_parts_old; ++i) {
/* Is the mass of this particle larger than the threshold? */
struct part *p = &parts[i];
/* Ignore inhibited particles */
if (part_is_inhibited(p, e)) continue;
const float gas_mass = hydro_get_mass(p);
const float h = p->h;
/* Found a particle to split */
if (gas_mass > mass_threshold) {
struct xpart *xp = &xparts[i];
struct gpart *gp = p->gpart;
/* Start by copying over the particles */
memcpy(&parts[k_parts], p, sizeof(struct part));
memcpy(&xparts[k_parts], xp, sizeof(struct xpart));
if (with_gravity) {
memcpy(&gparts[k_gparts], gp, sizeof(struct gpart));
}
/* Update the IDs */
parts[k_parts].id += 2;
/* Re-link everything */
if (with_gravity) {
parts[k_parts].gpart = &gparts[k_gparts];
gparts[k_gparts].id_or_neg_offset = -k_parts;
}
/* Displacement unit vector */
const double delta_x = random_unit_interval(p->id, e->ti_current,
(enum random_number_type)0);
const double delta_y = random_unit_interval(p->id, e->ti_current,
(enum random_number_type)1);
const double delta_z = random_unit_interval(p->id, e->ti_current,
(enum random_number_type)2);
/* Displace the old particle */
p->x[0] += delta_x * displacement_factor * h;
p->x[1] += delta_y * displacement_factor * h;
p->x[2] += delta_z * displacement_factor * h;
if (with_gravity) {
gp->x[0] += delta_x * displacement_factor * h;
gp->x[1] += delta_y * displacement_factor * h;
gp->x[2] += delta_z * displacement_factor * h;
}
/* Displace the new particle */
parts[k_parts].x[0] -= delta_x * displacement_factor * h;
parts[k_parts].x[1] -= delta_y * displacement_factor * h;
parts[k_parts].x[2] -= delta_z * displacement_factor * h;
if (with_gravity) {
gparts[k_gparts].x[0] -= delta_x * displacement_factor * h;
gparts[k_gparts].x[1] -= delta_y * displacement_factor * h;
gparts[k_gparts].x[2] -= delta_z * displacement_factor * h;
}
/* Divide the mass */
const double new_mass = gas_mass * 0.5;
hydro_set_mass(p, new_mass);
hydro_set_mass(&parts[k_parts], new_mass);
if (with_gravity) {
gp->mass = new_mass;
gparts[k_gparts].mass = new_mass;
}
/* Split the chemistry fields */
chemistry_split_part(p, particle_split_factor);
chemistry_split_part(&parts[k_parts], particle_split_factor);
/* Split the cooling fields */
cooling_split_part(p, xp, particle_split_factor);
cooling_split_part(&parts[k_parts], &xparts[k_parts],
particle_split_factor);
/* Split the star formation fields */
star_formation_split_part(p, xp, particle_split_factor);
star_formation_split_part(&parts[k_parts], &xparts[k_parts],
particle_split_factor);
/* Split the tracer fields */
tracers_split_part(p, xp, particle_split_factor);
tracers_split_part(&parts[k_parts], &xparts[k_parts],
particle_split_factor);
/* Mark the particles as not having been swallowed */
black_holes_mark_part_as_not_swallowed(&p->black_holes_data);
black_holes_mark_part_as_not_swallowed(&parts[k_parts].black_holes_data);
/* Move to the next free slot */
k_parts++;
if (with_gravity) k_gparts++;
}
}
/* Update the local counters */
s->nr_parts = k_parts;
s->nr_gparts = k_gparts;
#ifdef SWIFT_DEBUG_CHECKS
if (s->nr_parts != nr_parts_old + (particle_split_factor - 1) * counter) {
error("Incorrect number of particles created");
}
#endif
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
clocks_getunit());
}
......@@ -59,6 +59,8 @@ void hydro_props_init(struct hydro_props *p,
const struct unit_system *us,
struct swift_params *params) {
/* ------ Smoothing lengths parameters ---------- */
/* Kernel properties */
p->eta_neighbours = parser_get_param_float(params, "SPH:resolution_eta");
......@@ -98,6 +100,8 @@ void hydro_props_init(struct hydro_props *p,
if (p->max_smoothing_iterations <= 10)
error("The number of smoothing length iterations should be > 10");
/* ------ Neighbour number definition ------------ */
/* Non-conventional neighbour number definition */
p->use_mass_weighted_num_ngb =
parser_get_opt_param_int(params, "SPH:use_mass_weighted_num_ngb", 0);
......@@ -108,6 +112,8 @@ void hydro_props_init(struct hydro_props *p,
#endif
}
/* ------ Time integration parameters ------------ */
/* Time integration properties */
p->CFL_condition = parser_get_param_float(params, "SPH:CFL_condition");
const float max_volume_change = parser_get_opt_param_float(
......@@ -121,6 +127,8 @@ void hydro_props_init(struct hydro_props *p,
if (p->initial_temperature < 0.f)
error("ERROR: Initial temperature set to a negative value!!!");
/* ------ Temperature parameters ----------------- */
/* Minimal temperature */
p->minimal_temperature = parser_get_opt_param_float(
params, "SPH:minimal_temperature", hydro_props_default_min_temp);
......@@ -185,6 +193,17 @@ void hydro_props_init(struct hydro_props *p,
mean_molecular_weight = 4. / (1. + 3. * p->hydrogen_mass_fraction);
p->minimal_internal_energy = u_min / mean_molecular_weight;
/* ------ Particle splitting parameters ---------- */
/* Are we doing particle splitting? */
p->particle_splitting =
parser_get_opt_param_int(params, "SPH:particle_splitting", 0);
if (p->particle_splitting) {
p->particle_splitting_mass_threshold =
parser_get_param_float(params, "SPH:particle_splitting_mass_threshold");
}
}
/**
......@@ -238,6 +257,12 @@ void hydro_props_print(const struct hydro_props *p) {
if (p->minimal_temperature != hydro_props_default_min_temp)
message("Minimal gas temperature set to %f", p->minimal_temperature);
if (p->particle_splitting)
message("Splitting particles with mass > %e U_M",
p->particle_splitting_mass_threshold);
else
message("No particle splitting");
/* Print out the implementation-dependent viscosity parameters
* (see hydro/SCHEME/hydro_parameters.h for this implementation) */
viscosity_print(&(p->viscosity));
......
......@@ -47,6 +47,8 @@ struct unit_system;
*/
struct hydro_props {
/* ------ Smoothing lengths parameters ---------- */
/*! Resolution parameter */
float eta_neighbours;
......@@ -59,37 +61,43 @@ struct hydro_props {
/*! Tolerance on neighbour number (for info only)*/
float delta_neighbours;
/*! Maximal smoothing length */
/*! Maximal smoothing length (internal units) */
float h_max;
/*! Minimal smoothing length expressed as ratio to softening length */
float h_min_ratio;
/*! Minimal smoothing length */
/*! Minimal smoothing length (internal units) */
float h_min;
/*! Maximal number of iterations to converge h */
int max_smoothing_iterations;
/* ------ Neighbour number definition ------------ */
/*! Are we using the mass-weighted definition of neighbour number? */
int use_mass_weighted_num_ngb;
/* ------ Time integration parameters ------------ */
/*! Time integration properties */
float CFL_condition;
/*! Maximal change of h over one time-step */
float log_max_h_change;
/* ------ Temperature parameters ----------------- */
/*! Minimal temperature allowed */
float minimal_temperature;
/*! Minimal physical internal energy per unit mass */
/*! Minimal physical internal energy per unit mass (internal units) */
float minimal_internal_energy;
/*! Initial temperature */
float initial_temperature;
/*! Initial physical internal energy per unit mass */
/*! Initial physical internal energy per unit mass (internal units) */
float initial_internal_energy;
/*! Primordial hydrogen mass fraction for initial energy conversion */
......@@ -104,6 +112,16 @@ struct hydro_props {
/*! Mean molecular weight above hydrogen ionization temperature */
float mu_ionised;
/* ------ Particle splitting parameters ---------- */
/*! Is particle splitting activated? */
int particle_splitting;
/*! Mass above which particles get split (internal units) */
float particle_splitting_mass_threshold;
/* ------ Viscosity and diffusion ---------------- */
/*! Artificial viscosity parameters */
struct viscosity_global_data viscosity;
......
......@@ -722,4 +722,20 @@ star_formation_first_init_part(const struct phys_const* restrict phys_const,
const struct part* restrict p,
struct xpart* restrict xp) {}
/**
* @brief Split the star formation content of a particle into n pieces
*
* We only need to split the SFR if it is positive, i.e. it is not
* storing the redshift/time of last SF event.
*
* @param p The #part.
* @param xp The #xpart.
* @param n The number of pieces to split into.
*/
__attribute__((always_inline)) INLINE static void star_formation_split_part(
struct part* p, struct xpart* xp, const double n) {
if (xp->sf_data.SFR > 0.) xp->sf_data.SFR /= n;
}
#endif /* SWIFT_EAGLE_STAR_FORMATION_H */
......@@ -25,7 +25,8 @@
*/
struct star_formation_xpart_data {
/*! Star formation rate */
/*! Star formation rate (internal units) or (if negative) time/scale-factor of
* last SF episode */
float SFR;
};
......
......@@ -287,4 +287,17 @@ star_formation_first_init_part(const struct phys_const* restrict phys_const,
const struct part* restrict p,
struct xpart* restrict xp) {}
/**
* @brief Split the star formation content of a particle into n pieces
*
*
* @param p The #part.
* @param xp The #xpart.
* @param n The number of pieces to split into.
*/
__attribute__((always_inline)) INLINE static void star_formation_split_part(
struct part* p, struct xpart* xp, const double n) {
error("Loic: to be implemented");
}
#endif /* SWIFT_GEAR_STAR_FORMATION_H */
......@@ -227,4 +227,14 @@ star_formation_first_init_part(const struct phys_const* restrict phys_const,
const struct part* restrict p,
struct xpart* restrict xp) {}
/**
* @brief Split the star formation content of a particle into n pieces
*
* @param p The #part.
* @param xp The #xpart.
* @param n The number of pieces to split into.
*/
__attribute__((always_inline)) INLINE static void star_formation_split_part(
struct part* p, struct xpart* xp, const double n) {}
#endif /* SWIFT_NONE_STAR_FORMATION_H */
......@@ -158,4 +158,16 @@ static INLINE void tracers_after_black_holes_feedback(struct xpart *xp) {
xp->tracers_data.hit_by_AGN_feedback = 1;
}
/**
* @brief Split the tracer content of a particle into n pieces
*
* Nothing to do here.
*
* @param p The #part.
* @param xp The #xpart.
* @param n The number of pieces to split into.
*/
__attribute__((always_inline)) INLINE static void tracers_split_part(
struct part *p, struct xpart *xp, const double n) {}
#endif /* SWIFT_TRACERS_EAGLE_H */
......@@ -127,4 +127,15 @@ static INLINE void tracers_after_feedback(struct xpart *xp) {}
*/
static INLINE void tracers_after_black_holes_feedback(struct xpart *xp) {}
/**
* @brief Split the tracer content of a particle into n pieces
*
* Nothing to do here.
*
* @param p The #part.
* @param xp The #xpart.
* @param n The number of pieces to split into.
*/
__attribute__((always_inline)) INLINE static void tracers_split_part(
struct part *p, struct xpart *xp, const double n) {}
#endif /* SWIFT_TRACERS_NONE_H */
......@@ -54,6 +54,7 @@ threshold = 0.008
num_files = len(sys.argv) - 1
labels = [
["engine_split_gas_particles:", 1],
["Gpart assignment", 1],
["Mesh comunication", 1],
["Forward Fourier transform", 1],
......
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