Commit afbf0b3e authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'sink2' into 'master'

Implement sinks

See merge request !1069
parents 21144edc 345c50a5
......@@ -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
......@@ -105,6 +105,7 @@ int main(int argc, char *argv[]) {
struct space s;
struct spart *sparts = NULL;
struct bpart *bparts = NULL;
struct sink *sinks = NULL;
struct unit_system us;
struct los_props los_properties;
......@@ -171,6 +172,7 @@ int main(int argc, char *argv[]) {
int with_mpole_reconstruction = 0;
int with_structure_finding = 0;
int with_logger = 0;
int with_sink = 0;
int with_qla = 0;
int with_eagle = 0;
int with_gear = 0;
......@@ -221,6 +223,8 @@ int main(int argc, char *argv[]) {
OPT_BOOLEAN('S', "stars", &with_stars, "Run with stars.", NULL, 0, 0),
OPT_BOOLEAN('B', "black-holes", &with_black_holes,
"Run with black holes.", NULL, 0, 0),
OPT_BOOLEAN('k', "sinks", &with_sink, "Run with sink particles.", NULL, 0,
0),
OPT_BOOLEAN(
'u', "fof", &with_fof,
"Run Friends-of-Friends algorithm to perform black hole seeding.",
......@@ -439,6 +443,13 @@ int main(int argc, char *argv[]) {
}
#endif
#ifdef WITH_MPI
if (with_sink) {
printf("Error: sink particles are not available yet with MPI.\n");
return 1;
}
#endif
/* The CPU frequency is a long long, so we need to parse that ourselves. */
if (cpufreqarg != NULL) {
if (sscanf(cpufreqarg, "%llu", &cpufreq) != 1) {
......@@ -623,6 +634,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(bpart) is %4zi bytes.", sizeof(struct bpart));
message("sizeof(gpart) is %4zi bytes.", sizeof(struct gpart));
......@@ -1004,33 +1016,37 @@ int main(int argc, char *argv[]) {
fflush(stdout);
/* Get ready to read particles of all kinds */
size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0, Nspart = 0, Nbpart = 0;
size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0;
size_t Nspart = 0, Nbpart = 0, Nsink = 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_gravity, 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);
read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sinks, &sparts,
&bparts, &Ngas, &Ngpart, &Ngpart_background, &Nsink,
&Nspart, &Nbpart, &flag_entropy_ICs, with_hydro,
with_gravity, with_sink, 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);
#else
read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
&Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
&flag_entropy_ICs, with_hydro, with_gravity, 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);
read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sinks, &sparts,
&bparts, &Ngas, &Ngpart, &Ngpart_background, &Nsink, &Nspart,
&Nbpart, &flag_entropy_ICs, with_hydro, with_gravity,
with_sink, 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);
#endif
#else
read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
&Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
&flag_entropy_ICs, with_hydro, with_gravity, with_stars,
with_black_holes, with_cosmology, cleanup_h, cleanup_sqrt_a,
cosmo.h, cosmo.a, nr_threads, dry_run);
read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sinks, &sparts,
&bparts, &Ngas, &Ngpart, &Ngpart_background, &Nsink, &Nspart,
&Nbpart, &flag_entropy_ICs, with_hydro, with_gravity,
with_sink, with_stars, with_black_holes, with_cosmology,
cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads,
dry_run);
#endif
#endif
......@@ -1061,22 +1077,27 @@ int main(int argc, char *argv[]) {
for (size_t k = 0; k < Ngpart; ++k)
if (gparts[k].type == swift_type_gas) error("Linking problem");
}
if (!with_sink && !dry_run) {
for (size_t k = 0; k < Ngpart; ++k)
if (gparts[k].type == swift_type_sink) error("Linking problem");
}
/* Check that the other links are correctly set */
if (!dry_run)
part_verify_links(parts, gparts, sparts, bparts, Ngas, Ngpart, Nspart,
Nbpart, 1);
part_verify_links(parts, gparts, sinks, sparts, bparts, Ngas, Ngpart,
Nsink, Nspart, Nbpart, 1);
#endif
/* Get the total number of particles across all nodes. */
long long N_total[swift_type_count + 1] = {0};
long long Nbaryons = Ngas + Nspart + Nbpart;
long long Nbaryons = Ngas + Nspart + Nbpart + Nsink;
#if defined(WITH_MPI)
long long N_long[swift_type_count + 1] = {0};
N_long[swift_type_gas] = Ngas;
N_long[swift_type_dark_matter] =
with_gravity ? Ngpart - Ngpart_background - Nbaryons : 0;
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;
......@@ -1087,6 +1108,7 @@ int main(int argc, char *argv[]) {
N_total[swift_type_dark_matter] =
with_gravity ? Ngpart - Ngpart_background - Nbaryons : 0;
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;
......@@ -1094,17 +1116,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] >
N_total[swift_type_black_hole] + N_total[swift_type_sink] >
0) ||
(with_DM_particles && generate_gas_in_ics);
......@@ -1114,8 +1138,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_properties, parts, gparts,
sparts, bparts, Ngas, Ngpart, Nspart, Nbpart, periodic,
space_init(&s, params, &cosmo, dim, &hydro_properties, parts, gparts, sinks,
sparts, bparts, Ngas, Ngpart, Nsink, Nspart, Nbpart, periodic,
replicate, generate_gas_in_ics, with_hydro, with_self_gravity,
with_star_formation, with_DM_background_particles, talking,
dry_run, nr_nodes);
......@@ -1163,13 +1187,14 @@ int main(int argc, char *argv[]) {
space_check_cosmology(&s, &cosmo, myrank);
/* Also update the total counts (in case of changes due to replication) */
Nbaryons = s.nr_parts + s.nr_sparts + s.nr_bparts;
Nbaryons = s.nr_parts + s.nr_sparts + s.nr_bparts + s.nr_sinks;
#if defined(WITH_MPI)
N_long[swift_type_gas] = s.nr_parts;
N_long[swift_type_dark_matter] =
with_gravity ? s.nr_gparts - Ngpart_background - Nbaryons : 0;
N_long[swift_type_count] = s.nr_gparts;
N_long[swift_type_stars] = s.nr_sparts;
N_long[swift_type_sink] = s.nr_sinks;
N_long[swift_type_black_hole] = s.nr_bparts;
MPI_Allreduce(&N_long, &N_total, swift_type_count + 1, MPI_LONG_LONG_INT,
MPI_SUM, MPI_COMM_WORLD);
......@@ -1179,6 +1204,7 @@ int main(int argc, char *argv[]) {
with_gravity ? s.nr_gparts - Ngpart_background - Nbaryons : 0;
N_total[swift_type_count] = s.nr_gparts;
N_total[swift_type_stars] = s.nr_sparts;
N_total[swift_type_sink] = s.nr_sinks;
N_total[swift_type_black_hole] = s.nr_bparts;
#endif
......@@ -1191,6 +1217,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("%zi bparts in %i cells.", s.nr_bparts, s.tot_cells);
message("maximum depth is %d.", s.maxdepth);
......@@ -1247,12 +1274,13 @@ int main(int argc, char *argv[]) {
if (with_fof) engine_policies |= engine_policy_fof;
if (with_logger) engine_policies |= engine_policy_logger;
if (with_line_of_sight) engine_policies |= engine_policy_line_of_sight;
if (with_sink) engine_policies |= engine_policy_sink;
/* 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, &entropy_floor, &gravity_properties,
......@@ -1279,10 +1307,13 @@ 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 on %lld gas 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]);
"Running 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_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)...",
......
......@@ -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,