Commit 345c50a5 authored by Loic Hausammann's avatar Loic Hausammann Committed by Matthieu Schaller

Implement sinks

parent 21144edc
......@@ -784,6 +784,7 @@ INPUT += @top_srcdir@/src/stars/GEAR
INPUT += @top_srcdir@/src/feedback/EAGLE
INPUT += @top_srcdir@/src/feedback/GEAR
INPUT += @top_srcdir@/src/black_holes/EAGLE
INPUT += @top_srcdir@/src/sink/Default
INPUT += @top_srcdir@/logger
# This tag can be used to specify the character encoding of the source files
......
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/3e11-star-only-DM-halo-galaxy.hdf5
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.9891E43 # 10^10 solar masses
UnitLength_in_cgs: 3.08567758E21 # 1 kpc
UnitVelocity_in_cgs: 1E5 # km/s
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion).
max_physical_DM_softening: 0.7 # Physical softening length (in internal units).
max_physical_baryon_softening: 0.2 # Physical softening length (in internal units)
# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 1. # The end time of the simulation (in internal units).
dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: output # Common part of the name of output files
time_first: 0. # (Optional) Time of the first output if non-cosmological time-integration (in internal units)
delta_time: 0.01 # Time difference between consecutive outputs (in internal units)
Scheduler:
max_top_level_cells: 96
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-1 # Time between statistics output
time_first: 0. # (Optional) Time of the first stats output if non-cosmological time-integration (in internal units)
# Parameters related to the initial conditions
InitialConditions:
file_name: isolated_galaxy.hdf5 # The file to read
periodic: 0 # Are we running with periodic ICs?
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 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 numpy as np
from shutil import copyfile
# Add sink particles to the isolated galaxy
fileName = "3e11-star-only-DM-halo-galaxy.hdf5"
output = "isolated_galaxy.hdf5"
L = 13756 # Size of the box (internal units)
L_sink = 1000 # Size of the sink particle area (L_sink < L)
N = 1000 # Number of sink particles
max_vel = 100 # Maximal velocity of the sink particles (in km / s)
mass = 0.000142 # Mass of a sink particle (internal units)
min_id = 360001 # minimal ids of the other particles
# ---------------------------------------------------
copyfile(fileName, output)
# File
file = h5py.File(output, 'a')
pos = np.random.rand(N, 3) * L_sink + 0.5 * (L - L_sink)
vel = 2 * (np .random.rand(N, 3) - 0.5) * max_vel
m = mass * np.ones(N)
# Generate extra arrays
ids = np.linspace(min_id, min_id + N, N)
# --------------------------------------------------
# Header
grp = file["Header"]
numpart = grp.attrs["NumPart_Total"][:]
numpart[3] = N
grp.attrs["NumPart_Total"] = numpart
grp.attrs["NumPart_ThisFile"] = numpart
# Units
grp = file["Units"]
uv = grp.attrs["Unit length in cgs (U_L)"]
uv /= grp.attrs["Unit time in cgs (U_t)"]
vel *= 1e5 / uv
# Particle group
grp = file.create_group("/PartType3")
grp.create_dataset('Coordinates', data=pos, dtype='d')
grp.create_dataset('Velocities', data=vel, dtype='f')
grp.create_dataset('Masses', data=m, dtype='f')
grp.create_dataset('ParticleIDs', data=ids, dtype='L')
file.close()
#!/bin/bash
if [ ! -e 3e11-star-only-DM-halo-galaxy.hdf5 ]
then
echo "Fetching initial conditons for the isolated galaxy with an external potential ..."
./getIC.sh
fi
echo "Generate initial conditions"
python3 makeIC.py
../../swift --sinks --external-gravity --self-gravity --threads=16 isolated_galaxy.yml 2>&1 | tee output.log
This diff is collapsed.
......@@ -95,6 +95,7 @@ int main(int argc, char *argv[]) {
struct phys_const prog_const;
struct space s;
struct spart *sparts = NULL;
struct sink *sinks = NULL;
struct bpart *bparts = NULL;
struct unit_system us;
......@@ -141,6 +142,7 @@ int main(int argc, char *argv[]) {
int dump_threadpool = 0;
int with_fp_exceptions = 0;
int with_cosmology = 0;
int with_sinks = 0;
int with_stars = 0;
int with_black_holes = 0;
int with_hydro = 0;
......@@ -164,6 +166,8 @@ int main(int argc, char *argv[]) {
"Run with cosmological information.", NULL, 0, 0),
OPT_BOOLEAN(0, "hydro", &with_hydro, "Read gas particles from the ICs.",
NULL, 0, 0),
OPT_BOOLEAN(0, "sinks", &with_sinks, "Read sinks from the ICs.", NULL, 0,
0),
OPT_BOOLEAN(0, "stars", &with_stars, "Read stars from the ICs.", NULL, 0,
0),
OPT_BOOLEAN(0, "black-holes", &with_black_holes,
......@@ -317,6 +321,7 @@ int main(int argc, char *argv[]) {
if (myrank == 0) {
message("sizeof(part) is %4zi bytes.", sizeof(struct part));
message("sizeof(xpart) is %4zi bytes.", sizeof(struct xpart));
message("sizeof(sink) is %4zi bytes.", sizeof(struct sink));
message("sizeof(spart) is %4zi bytes.", sizeof(struct spart));
message("sizeof(gpart) is %4zi bytes.", sizeof(struct gpart));
message("sizeof(multipole) is %4zi bytes.", sizeof(struct multipole));
......@@ -447,34 +452,37 @@ int main(int argc, char *argv[]) {
/* Get ready to read particles of all kinds */
int flag_entropy_ICs = 0;
size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0, Nspart = 0, Nbpart = 0;
size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0;
size_t Nsink = 0, Nspart = 0, Nbpart = 0;
double dim[3] = {0., 0., 0.};
if (myrank == 0) clocks_gettime(&tic);
#if defined(HAVE_HDF5)
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
&Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
&flag_entropy_ICs, with_hydro,
/*with_grav=*/1, with_stars, with_black_holes,
read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sinks, &sparts,
&bparts, &Ngas, &Ngpart, &Ngpart_background, &Nsink, &Nspart,
&Nbpart, &flag_entropy_ICs, with_hydro,
/*with_grav=*/1, with_sinks, with_stars, with_black_holes,
with_cosmology, cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a,
myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
/*dry_run=*/0);
#else
read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas,
&Ngpart, &Ngpart_background, &Nspart, &Nbpart,
&flag_entropy_ICs, with_hydro,
/*with_grav=*/1, with_stars, with_black_holes, with_cosmology,
cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank, nr_nodes,
MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, /*dry_run=*/0);
read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sinks, &sparts,
&bparts, &Ngas, &Ngpart, &Ngpart_background, &Nsink, &Nspart,
&Nbpart, &flag_entropy_ICs, with_hydro,
/*with_grav=*/1, with_sinks, with_stars, with_black_holes,
with_cosmology, cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a,
myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
/*dry_run=*/0);
#endif
#else
read_ic_single(
ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas, &Ngpart,
&Ngpart_background, &Nspart, &Nbpart, &flag_entropy_ICs, with_hydro,
/*with_grav=*/1, with_stars, with_black_holes, with_cosmology, cleanup_h,
cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads, /*dry_run=*/0);
read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sinks, &sparts,
&bparts, &Ngas, &Ngpart, &Ngpart_background, &Nsink, &Nspart,
&Nbpart, &flag_entropy_ICs, with_hydro,
/*with_grav=*/1, with_sinks, with_stars, with_black_holes,
with_cosmology, cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a,
nr_threads, /*dry_run=*/0);
#endif
#endif
if (myrank == 0) {
......@@ -498,6 +506,7 @@ int main(int argc, char *argv[]) {
N_long[swift_type_gas] = Ngas;
N_long[swift_type_dark_matter] = Ngpart - Ngpart_background - Nbaryons;
N_long[swift_type_dark_matter_background] = Ngpart_background;
N_long[swift_type_sink] = Nsink;
N_long[swift_type_stars] = Nspart;
N_long[swift_type_black_hole] = Nbpart;
N_long[swift_type_count] = Ngpart;
......@@ -507,6 +516,7 @@ int main(int argc, char *argv[]) {
N_total[swift_type_gas] = Ngas;
N_total[swift_type_dark_matter] = Ngpart - Ngpart_background - Nbaryons;
N_total[swift_type_dark_matter_background] = Ngpart_background;
N_total[swift_type_sink] = Nsink;
N_total[swift_type_stars] = Nspart;
N_total[swift_type_black_hole] = Nbpart;
N_total[swift_type_count] = Ngpart;
......@@ -514,17 +524,19 @@ int main(int argc, char *argv[]) {
if (myrank == 0)
message(
"Read %lld gas particles, %lld stars particles, %lld black hole "
"Read %lld gas particles, %lld sink particles, %lld stars particles, "
"%lld black hole "
"particles, %lld DM particles and %lld DM background particles from "
"the ICs.",
N_total[swift_type_gas], N_total[swift_type_stars],
N_total[swift_type_black_hole], N_total[swift_type_dark_matter],
N_total[swift_type_gas], N_total[swift_type_sink],
N_total[swift_type_stars], N_total[swift_type_black_hole],
N_total[swift_type_dark_matter],
N_total[swift_type_dark_matter_background]);
const int with_DM_particles = N_total[swift_type_dark_matter] > 0;
const int with_baryon_particles =
(N_total[swift_type_gas] + N_total[swift_type_stars] +
N_total[swift_type_black_hole]) > 0;
(N_total[swift_type_gas] + N_total[swift_type_sink] +
N_total[swift_type_stars] + N_total[swift_type_black_hole]) > 0;
/* Do we have background DM particles? */
const int with_DM_background_particles =
......@@ -533,7 +545,8 @@ int main(int argc, char *argv[]) {
/* Initialize the space with these data. */
if (myrank == 0) clocks_gettime(&tic);
space_init(&s, params, &cosmo, dim, /*hydro_props=*/NULL, parts, gparts,
sparts, bparts, Ngas, Ngpart, Nspart, Nbpart, periodic, replicate,
sinks, sparts, bparts, Ngas, Ngpart, Nsink, Nspart, Nbpart,
periodic, replicate,
/*generate_gas_in_ics=*/0, /*hydro=*/N_total[0] > 0, /*gravity=*/1,
/*with_star_formation=*/0, with_DM_background_particles, talking,
/*dry_run=*/0, nr_nodes);
......@@ -570,6 +583,7 @@ int main(int argc, char *argv[]) {
N_long[swift_type_gas] = s.nr_parts;
N_long[swift_type_dark_matter] = s.nr_gparts - Ngpart_background - Nbaryons;
N_long[swift_type_count] = s.nr_gparts;
N_long[swift_type_sink] = s.nr_sinks;
N_long[swift_type_stars] = s.nr_sparts;
N_long[swift_type_black_hole] = s.nr_bparts;
MPI_Allreduce(&N_long, &N_total, swift_type_count + 1, MPI_LONG_LONG_INT,
......@@ -578,6 +592,7 @@ int main(int argc, char *argv[]) {
N_total[swift_type_gas] = s.nr_parts;
N_total[swift_type_dark_matter] = s.nr_gparts - Ngpart_background - Nbaryons;
N_total[swift_type_count] = s.nr_gparts;
N_total[swift_type_sink] = s.nr_sinks;
N_total[swift_type_stars] = s.nr_sparts;
N_total[swift_type_black_hole] = s.nr_bparts;
#endif
......@@ -591,6 +606,7 @@ int main(int argc, char *argv[]) {
s.cdim[1], s.cdim[2]);
message("%zi parts in %i cells.", s.nr_parts, s.tot_cells);
message("%zi gparts in %i cells.", s.nr_gparts, s.tot_cells);
message("%zi sinks in %i cells.", s.nr_sinks, s.tot_cells);
message("%zi sparts in %i cells.", s.nr_sparts, s.tot_cells);
message("maximum depth is %d.", s.maxdepth);
fflush(stdout);
......@@ -605,8 +621,8 @@ int main(int argc, char *argv[]) {
/* Initialize the engine with the space and policies. */
if (myrank == 0) clocks_gettime(&tic);
engine_init(&e, &s, params, output_options, N_total[swift_type_gas],
N_total[swift_type_count], N_total[swift_type_stars],
N_total[swift_type_black_hole],
N_total[swift_type_count], N_total[swift_type_sink],
N_total[swift_type_stars], N_total[swift_type_black_hole],
N_total[swift_type_dark_matter_background], engine_policies,
talking, &reparttype, &us, &prog_const, &cosmo,
/*hydro_properties=*/NULL, /*entropy_floor=*/NULL,
......@@ -630,10 +646,12 @@ int main(int argc, char *argv[]) {
const long long N_DM = N_total[swift_type_dark_matter] +
N_total[swift_type_dark_matter_background];
message(
"Running FOF on %lld gas particles, %lld stars particles %lld black "
"Running FOF on %lld gas particles, %lld sink particles, %lld stars "
"particles %lld black "
"hole particles and %lld DM particles (%lld gravity particles)",
N_total[swift_type_gas], N_total[swift_type_stars],
N_total[swift_type_black_hole], N_DM, N_total[swift_type_count]);
N_total[swift_type_gas], N_total[swift_type_sink],
N_total[swift_type_stars], N_total[swift_type_black_hole], N_DM,
N_total[swift_type_count]);
message(
"from t=%.3e until t=%.3e with %d ranks, %d threads / rank and %d "
"task queues / rank (dt_min=%.3e, dt_max=%.3e)...",
......
......@@ -1648,6 +1648,8 @@ void cell_bunlocktree(struct cell *c) {
* @param bparts_offset Offset of the cell bparts array relative to the
* space's bparts array, i.e. c->black_holes.parts -
* s->black_holes.parts.
* @param sinks_offset Offset of the cell sink array relative to the
* space's sink array, i.e. c->sinks.parts - s->sinks.parts.
* @param buff A buffer with at least max(c->hydro.count, c->grav.count)
* entries, used for sorting indices.
* @param sbuff A buffer with at least max(c->stars.count, c->grav.count)
......@@ -1656,18 +1658,23 @@ void cell_bunlocktree(struct cell *c) {
* entries, used for sorting indices for the sparts.
* @param gbuff A buffer with at least max(c->hydro.count, c->grav.count)
* entries, used for sorting indices for the gparts.
* @param sinkbuff A buffer with at least max(c->sinks.count, c->grav.count)
* entries, used for sorting indices for the sinks.
*/
void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
ptrdiff_t bparts_offset, struct cell_buff *buff,
struct cell_buff *sbuff, struct cell_buff *bbuff,
struct cell_buff *gbuff) {
ptrdiff_t bparts_offset, ptrdiff_t sinks_offset,
struct cell_buff *buff, struct cell_buff *sbuff,
struct cell_buff *bbuff, struct cell_buff *gbuff,
struct cell_buff *sinkbuff) {
const int count = c->hydro.count, gcount = c->grav.count,
scount = c->stars.count, bcount = c->black_holes.count;
scount = c->stars.count, bcount = c->black_holes.count,
sink_count = c->sinks.count;
struct part *parts = c->hydro.parts;
struct xpart *xparts = c->hydro.xparts;
struct gpart *gparts = c->grav.parts;
struct spart *sparts = c->stars.parts;
struct bpart *bparts = c->black_holes.parts;
struct sink *sinks = c->sinks.parts;
const double pivot[3] = {c->loc[0] + c->width[0] / 2,
c->loc[1] + c->width[1] / 2,
c->loc[2] + c->width[2] / 2};
......@@ -1696,6 +1703,11 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
bbuff[k].x[2] != bparts[k].x[2])
error("Inconsistent bbuff contents.");
}
for (int k = 0; k < sink_count; k++) {
if (sinkbuff[k].x[0] != sinks[k].x[0] ||
sinkbuff[k].x[1] != sinks[k].x[1] || sinkbuff[k].x[2] != sinks[k].x[2])
error("Inconsistent sinkbuff contents.");
}
#endif /* SWIFT_DEBUG_CHECKS */
/* Fill the buffer with the indices. */
......@@ -1925,6 +1937,60 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
c->progeny[k]->black_holes.parts = &c->black_holes.parts[bucket_offset[k]];
}
/* Now do the same song and dance for the sinks. */
for (int k = 0; k < 8; k++) bucket_count[k] = 0;
/* Fill the buffer with the indices. */
for (int k = 0; k < sink_count; k++) {
const int bid = (sinkbuff[k].x[0] > pivot[0]) * 4 +
(sinkbuff[k].x[1] > pivot[1]) * 2 +
(sinkbuff[k].x[2] > pivot[2]);
bucket_count[bid]++;
sinkbuff[k].ind = bid;
}
/* Set the buffer offsets. */
bucket_offset[0] = 0;
for (int k = 1; k <= 8; k++) {
bucket_offset[k] = bucket_offset[k - 1] + bucket_count[k - 1];
bucket_count[k - 1] = 0;
}
/* Run through the buckets, and swap particles to their correct spot. */
for (int bucket = 0; bucket < 8; bucket++) {
for (int k = bucket_offset[bucket] + bucket_count[bucket];
k < bucket_offset[bucket + 1]; k++) {
int bid = sinkbuff[k].ind;
if (bid != bucket) {
struct sink sink = sinks[k];
struct cell_buff temp_buff = sinkbuff[k];
while (bid != bucket) {
int j = bucket_offset[bid] + bucket_count[bid]++;
while (sinkbuff[j].ind == bid) {
j++;
bucket_count[bid]++;
}
memswap(&sinks[j], &sink, sizeof(struct sink));
memswap(&sinkbuff[j], &temp_buff, sizeof(struct cell_buff));
if (sinks[j].gpart)
sinks[j].gpart->id_or_neg_offset = -(j + sinks_offset);
bid = temp_buff.ind;
}
sinks[k] = sink;
sinkbuff[k] = temp_buff;
if (sinks[k].gpart)
sinks[k].gpart->id_or_neg_offset = -(k + sinks_offset);
}
bucket_count[bid]++;
}
}
/* Store the counts and offsets. */
for (int k = 0; k < 8; k++) {
c->progeny[k]->sinks.count = bucket_count[k];
c->progeny[k]->sinks.parts = &c->sinks.parts[bucket_offset[k]];
}
/* Finally, do the same song and dance for the gparts. */
for (int k = 0; k < 8; k++) bucket_count[k] = 0;
......@@ -1965,6 +2031,9 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
} else if (gparts[j].type == swift_type_stars) {
sparts[-gparts[j].id_or_neg_offset - sparts_offset].gpart =
&gparts[j];
} else if (gparts[j].type == swift_type_sink) {
sinks[-gparts[j].id_or_neg_offset - sinks_offset].gpart =
&gparts[j];
} else if (gparts[j].type == swift_type_black_hole) {
bparts[-gparts[j].id_or_neg_offset - bparts_offset].gpart =
&gparts[j];
......@@ -1978,6 +2047,8 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
} else if (gparts[k].type == swift_type_stars) {
sparts[-gparts[k].id_or_neg_offset - sparts_offset].gpart =
&gparts[k];
} else if (gparts[k].type == swift_type_sink) {
sinks[-gparts[k].id_or_neg_offset - sinks_offset].gpart = &gparts[k];
} else if (gparts[k].type == swift_type_black_hole) {
bparts[-gparts[k].id_or_neg_offset - bparts_offset].gpart =
&gparts[k];
......
......@@ -752,6 +752,19 @@ struct cell {
} black_holes;
/*! Sink particles variables */
struct {
/*! Pointer to the #sink data. */
struct sink *parts;
/*! Nr of #sink in this cell. */
int count;
/*! Nr of #sink this cell can hold after addition of new one. */
int count_total;
} sinks;
#ifdef WITH_MPI
/*! MPI variables */
struct {
......@@ -843,9 +856,10 @@ struct cell {
/* Function prototypes. */
void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
ptrdiff_t bparts_offset, struct cell_buff *buff,
struct cell_buff *sbuff, struct cell_buff *bbuff,
struct cell_buff *gbuff);
ptrdiff_t bparts_offset, ptrdiff_t sinks_offset,
struct cell_buff *buff, struct cell_buff *sbuff,
struct cell_buff *bbuff, struct cell_buff *gbuff,
struct cell_buff *sinkbuff);
void cell_sanitize(struct cell *c, int treated);
int cell_locktree(struct cell *c);
void cell_unlocktree(struct cell *c);
......
This diff is collapsed.
......@@ -40,6 +40,7 @@ struct velociraptor_gpart_data;
struct spart;
struct bpart;
struct xpart;
struct sink;
struct io_props;
struct engine;
struct threadpool;
......@@ -134,6 +135,10 @@ void io_collect_parts_to_write(const struct part* restrict parts,
struct xpart* restrict xparts_written,
const size_t Nparts,
const size_t Nparts_written);
void io_collect_sinks_to_write(const struct sink* restrict sinks,
struct sink* restrict sinks_written,
const size_t Nsinks,
const size_t Nsinks_written);
void io_collect_sparts_to_write(const struct spart* restrict sparts,
struct spart* restrict sparts_written,
const size_t Nsparts,
......@@ -167,6 +172,9 @@ void io_duplicate_stars_gparts(struct threadpool* tp,
struct spart* const sparts,
struct gpart* const gparts, size_t Nstars,
size_t Ndm);
void io_duplicate_sinks_gparts(struct threadpool* tp, struct sink* const sinks,
struct gpart* const gparts, size_t Nsinks,
size_t Ndm);
void io_duplicate_black_holes_gparts(struct threadpool* tp,
struct bpart* const bparts,
struct gpart* const gparts, size_t Nstars,
......
......@@ -54,6 +54,7 @@
#include "output_options.h"
#include "part.h"
#include "part_type.h"
#include "sink_io.h"
#include "star_formation_io.h"
#include "stars_io.h"
#include "tools.h"
......@@ -248,6 +249,7 @@ void write_output_distributed(struct engine* e,
const struct part* parts = e->s->parts;
const struct xpart* xparts = e->s->xparts;
const struct gpart* gparts = e->s->gparts;
const struct sink* sinks = e->s->sinks;
const struct spart* sparts = e->s->sparts;
const struct bpart* bparts = e->s->bparts;
struct output_options* output_options = e->output_options;
......@@ -267,6 +269,7 @@ void write_output_distributed(struct engine* e,
/* Number of particles currently in the arrays */
const size_t Ntot = e->s->nr_gparts;
const size_t Ngas = e->s->nr_parts;
const size_t Nsinks = e->s->nr_sinks;
const size_t Nstars = e->s->nr_sparts;
const size_t Nblackholes = e->s->nr_bparts;
......@@ -281,12 +284,14 @@ void write_output_distributed(struct engine* e,
e->s->nr_gparts - e->s->nr_inhibited_gparts - e->s->nr_extra_gparts;
const size_t Ngas_written =
e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
const size_t Nsinks_written =
e->s->nr_sinks - e->s->nr_inhibited_sinks - e->s->nr_extra_sinks;
const size_t Nstars_written =
e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
const size_t Nblackholes_written =
e->s->nr_bparts - e->s->nr_inhibited_bparts - e->s->nr_extra_bparts;
const size_t Nbaryons_written =
Ngas_written + Nstars_written + Nblackholes_written;
Ngas_written + Nstars_written + Nblackholes_written + Nsinks_written;
const size_t Ndm_written =
Ntot_written > 0 ? Ntot_written - Nbaryons_written - Ndm_background : 0;
......@@ -331,7 +336,7 @@ void write_output_distributed(struct engine* e,
/* Compute offset in the file and total number of particles */
const long long N[swift_type_count] = {Ngas_written, Ndm_written,
Ndm_background, 0,
Ndm_background, Nsinks_written,
Nstars_written, Nblackholes_written};
/* Gather the total number of particles to write */
......@@ -471,6 +476,7 @@ void write_output_distributed(struct engine* e,
struct xpart* xparts_written = NULL;
struct gpart* gparts_written = NULL;
struct velociraptor_gpart_data* gpart_group_data_written = NULL;
struct sink* sinks_written = NULL;
struct spart* sparts_written = NULL;
struct bpart* bparts_written = NULL;
......@@ -632,6 +638,34 @@ void write_output_distributed(struct engine* e,
}
} break;
case swift_type_sink: {
if (Nsinks == Nsinks_written) {
/* No inhibted particles: easy case */
Nparticles = Nsinks;
sink_write_particles(sinks, list, &num_fields, with_cosmology);
} else {
/* Ok, we need to fish out the particles we want */
Nparticles = Nsinks_written;
/* Allocate temporary arrays */
if (swift_memalign("sinks_written", (void**)&sinks_written,
sink_align,
Nsinks_written * sizeof(struct sink)) != 0)
error("Error while allocating temporary memory for sinks");
/* Collect the particles we want to write */
io_collect_sinks_to_write(sinks, sinks_written, Nsinks,
Nsinks_written);
/* Select the fields to write */
sink_write_particles(sinks_written, list, &num_fields,
with_cosmology);
}
} break;
case swift_type_stars: {
if (Nstars == Nstars_written) {
......@@ -759,6 +793,7 @@ void write_output_distributed(struct engine* e,
if (gparts_written) swift_free("gparts_written", gparts_written);
if (gpart_group_data_written)
swift_free("gpart_group_written", gpart_group_data_written);
if (sinks_written) swift_free("sinks_written", sinks_written);
if (sparts_written) swift_free("sparts_written", sparts_written);
if (bparts_written) swift_free("bparts_written", bparts_written);
......
......@@ -128,7 +128,8 @@ const char *engine_policy_names[] = {"none",
"time-step limiter",
"time-step sync",
"logger",
"line of sight"};
"line of sight",
"sink"};
/** The rank of the engine as a global variable (for messages). */
int engine_rank;
......@@ -1513,6 +1514,7 @@ void engine_print_task_counts(const struct engine *e) {
fflush(stdout);
message("nr_parts = %zu.", e->s->nr_parts);
message("nr_gparts = %zu.", e->s->nr_gparts);
message("nr_sink = %zu.", e->s->nr_sinks);
message("nr_sparts = %zu.", e->s->nr_sparts);
message("nr_bparts = %zu.", e->s->nr_bparts);
......@@ -1770,9 +1772,10 @@ void engine_rebuild(struct engine *e, const int repartitioned,
engine_recompute_displacement_constraint(e);
#ifdef SWIFT_DEBUG_CHECKS
part_verify_links(e->s->parts, e->s->gparts, e->s->sparts, e->s->bparts,
e->s->nr_parts, e->s->nr_gparts, e->s->nr_sparts,
e->s->nr_bparts, e->verbose);
part_verify_links(e->s->parts, e->s->gparts, e->s->sinks, e->s->sparts,
e->s->bparts, e->s->nr_parts, e->s->nr_gparts,
e->s->nr_sinks, e->s->nr_sparts, e->s->nr_bparts,
e->verbose);
#endif
/* Initial cleaning up session ? */
......@@ -2160,6 +2163,7 @@ void engine_first_init_particles(struct engine *e) {
space_first_init_gparts(e->s, e->verbose);
space_first_init_sparts(e->s, e->verbose);
space_first_init_bparts(e->s, e->verbose);
space_first_init_sinks(e->s, e->verbose);
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
......@@ -2213,6 +2217,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
space_init_gparts(s, e->verbose);
space_init_sparts(s, e->verbose);
space_init_bparts(s, e->verbose);
space_init_sinks(s, e->verbose);
/* Update the cooling function */
if ((e->policy & engine_policy_cooling) ||
......@@ -2284,6 +2289,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
space_init_gparts(e->s, e->verbose);
space_init_sparts(e->s, e->verbose);
space_init_bparts(e->s, e->verbose);
space_init_sinks(e->s, e->verbose);
/* Print the number of active tasks ? */
if (e->verbose) engine_print_task_counts(e);
......@@ -2419,9 +2425,10 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
#ifdef SWIFT_DEBUG_CHECKS
space_check_timesteps(e->s);
part_verify_links(e->s->parts, e->s->gparts, e->s->sparts, e->s->bparts,
e->s->nr_parts, e->s->nr_gparts, e->s->nr_sparts,
e->s->nr_bparts, e->verbose);
part_verify_links(e->s->parts, e->s->gparts, e->s->sinks, e->s->sparts,
e->s->bparts, e->s->nr_parts, e->s->nr_gparts,
e->s->nr_sinks, e->s->nr_sparts, e->s->nr_bparts,
e->verbose);
#endif
/* Ready to go */
......@@ -3483,14 +3490,15 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
/* Re-link everything to the gparts. */
if (s->nr_gparts > 0)
part_relink_all_parts_to_gparts(s->gparts, s->nr_gparts, s->parts,
part_relink_all_parts_to_gparts(s->gparts, s->nr_gparts, s->parts, s->sinks,
s<