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

Can now read and write star particles in single-node mode.

parent 982a3664
......@@ -762,6 +762,7 @@ WARN_LOGFILE =
INPUT = @top_srcdir@ @top_srcdir@/src @top_srcdir@/tests @top_srcdir@/examples
INPUT += @top_srcdir@/src/hydro/Minimal
INPUT += @top_srcdir@/src/gravity/Default
INPUT += @top_srcdir@/src/stars/Default
INPUT += @top_srcdir@/src/riemann
INPUT += @top_srcdir@/src/potential/point_mass
INPUT += @top_srcdir@/src/cooling/const_du
......
......@@ -303,6 +303,7 @@ int main(int argc, char *argv[]) {
if (myrank == 0) {
message("sizeof(struct part) is %4zi bytes.", sizeof(struct part));
message("sizeof(struct xpart) is %4zi bytes.", sizeof(struct xpart));
message("sizeof(struct spart) is %4zi bytes.", sizeof(struct spart));
message("sizeof(struct gpart) is %4zi bytes.", sizeof(struct gpart));
message("sizeof(struct task) is %4zi bytes.", sizeof(struct task));
message("sizeof(struct cell) is %4zi bytes.", sizeof(struct cell));
......@@ -368,7 +369,8 @@ int main(int argc, char *argv[]) {
struct part *parts = NULL;
struct gpart *gparts = NULL;
size_t Ngas = 0, Ngpart = 0;
struct spart *sparts = NULL;
size_t Ngas = 0, Ngpart = 0, Nspart = 0;
double dim[3] = {0., 0., 0.};
int periodic = 0;
int flag_entropy_ICs = 0;
......@@ -384,8 +386,8 @@ int main(int argc, char *argv[]) {
MPI_INFO_NULL, dry_run);
#endif
#else
read_ic_single(ICfileName, &us, dim, &parts, &gparts, &Ngas, &Ngpart,
&periodic, &flag_entropy_ICs, dry_run);
read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas, &Ngpart,
&Nspart, &periodic, &flag_entropy_ICs, dry_run);
#endif
if (myrank == 0) {
clocks_gettime(&toc);
......@@ -400,6 +402,7 @@ int main(int argc, char *argv[]) {
free(gparts);
gparts = NULL;
for (size_t k = 0; k < Ngas; ++k) parts[k].gpart = NULL;
for (size_t k = 0; k < Nspart; ++k) sparts[k].gpart = NULL;
Ngpart = 0;
}
if (!with_hydro) {
......@@ -411,23 +414,26 @@ int main(int argc, char *argv[]) {
}
/* Get the total number of particles across all nodes. */
long long N_total[2] = {0, 0};
long long N_total[3] = {0, 0, 0};
#if defined(WITH_MPI)
long long N_long[2] = {Ngas, Ngpart};
MPI_Reduce(&N_long, &N_total, 2, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
#else
N_total[0] = Ngas;
N_total[1] = Ngpart;
N_total[2] = Nspart;
#endif
if (myrank == 0)
message("Read %lld gas particles and %lld gparts from the ICs.", N_total[0],
N_total[1]);
message(
"Read %lld gas particles, %lld star particles and %lld gparts from the "
"ICs.",
N_total[0], N_total[2], N_total[1]);
/* Initialize the space with these data. */
if (myrank == 0) clocks_gettime(&tic);
struct space s;
space_init(&s, params, dim, parts, gparts, Ngas, Ngpart, periodic,
with_self_gravity, talking, dry_run);
space_init(&s, params, dim, parts, gparts, sparts, Ngas, Ngpart, Nspart,
periodic, with_self_gravity, talking, dry_run);
if (myrank == 0) {
clocks_gettime(&toc);
message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
......@@ -529,6 +535,8 @@ int main(int argc, char *argv[]) {
return 0;
}
engine_dump_snapshot(&e);
#ifdef WITH_MPI
/* Split the space. */
engine_split(&e, &initial_partition);
......
......@@ -597,7 +597,7 @@ void prepare_dm_gparts(struct gpart* const gparts, size_t Ndm) {
*
* @param parts The array of #part freshly read in.
* @param gparts The array of #gpart freshly read in with all the DM particles
*at the start
* at the start
* @param Ngas The number of gas particles read in.
* @param Ndm The number of DM particles read in.
*/
......@@ -624,6 +624,41 @@ void duplicate_hydro_gparts(struct part* const parts,
}
}
/**
* @brief Copy every #spart into the corresponding #gpart and link them.
*
* This function assumes that the DM particles and gas particles are all at
* the start of the gparts array and adds the star particles afterwards
*
* @param sparts The array of #spart freshly read in.
* @param gparts The array of #gpart freshly read in with all the DM and gas
* particles at the start.
* @param Nstars The number of stars particles read in.
* @param Ndm The number of DM and gas particles read in.
*/
void duplicate_star_gparts(struct spart* const sparts,
struct gpart* const gparts, size_t Nstars,
size_t Ndm) {
for (size_t i = 0; i < Nstars; ++i) {
/* Duplicate the crucial information */
gparts[i + Ndm].x[0] = sparts[i].x[0];
gparts[i + Ndm].x[1] = sparts[i].x[1];
gparts[i + Ndm].x[2] = sparts[i].x[2];
gparts[i + Ndm].v_full[0] = sparts[i].v[0];
gparts[i + Ndm].v_full[1] = sparts[i].v[1];
gparts[i + Ndm].v_full[2] = sparts[i].v[2];
gparts[i + Ndm].mass = sparts[i].mass;
/* Link the particles */
gparts[i + Ndm].id_or_neg_offset = -i;
sparts[i].gpart = &gparts[i + Ndm];
}
}
/**
* @brief Copy every DM #gpart into the dmparts array.
*
......
......@@ -75,6 +75,9 @@ void prepare_dm_gparts(struct gpart* const gparts, size_t Ndm);
void duplicate_hydro_gparts(struct part* const parts,
struct gpart* const gparts, size_t Ngas,
size_t Ndm);
void duplicate_star_gparts(struct spart* const sparts,
struct gpart* const gparts, size_t Nstars,
size_t Ndm);
void readAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data);
......
......@@ -36,6 +36,7 @@
/* Some constants. */
#define part_align 128
#define xpart_align 128
#define spart_align 128
#define gpart_align 128
/* Import the right hydro particle definition */
......@@ -62,6 +63,9 @@
/* Import the right gravity particle definition */
#include "./gravity/Default/gravity_part.h"
/* Import the right star particle definition */
#include "./stars/Default/star_part.h"
void part_relink_gparts(struct part *parts, size_t N, ptrdiff_t offset);
void part_relink_parts(struct gpart *gparts, size_t N, struct part *parts);
#ifdef WITH_MPI
......
......@@ -45,6 +45,7 @@
#include "io_properties.h"
#include "kernel_hydro.h"
#include "part.h"
#include "stars_io.h"
#include "units.h"
/*-----------------------------------------------------------------------------
......@@ -312,10 +313,12 @@ void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
* @param fileName The file to read.
* @param internal_units The system units used internally
* @param dim (output) The dimension of the volume.
* @param parts (output) Array of Gas particles.
* @param parts (output) Array of #part particles.
* @param gparts (output) Array of #gpart particles.
* @param sparts (output) Array of #spart particles.
* @param Ngas (output) number of Gas particles read.
* @param Ngparts (output) The number of #gpart read.
* @param Nstars (output) The number of #spart read.
* @param periodic (output) 1 if the volume is periodic, 0 if not.
* @param flag_entropy (output) 1 if the ICs contained Entropy in the
* InternalEnergy
......@@ -332,8 +335,9 @@ void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
*/
void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
size_t* Ngas, size_t* Ngparts, int* periodic,
int* flag_entropy, int dry_run) {
struct spart** sparts, size_t* Ngas, size_t* Ngparts,
size_t* Nstars, int* periodic, int* flag_entropy,
int dry_run) {
hid_t h_file = 0, h_grp = 0;
/* GADGET has only cubic boxes (in cosmological mode) */
......@@ -439,15 +443,22 @@ void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
units_conversion_factor(ic_units, internal_units, UNIT_CONV_LENGTH);
/* Allocate memory to store SPH particles */
*Ngas = N[0];
*Ngas = N[GAS];
if (posix_memalign((void*)parts, part_align, *Ngas * sizeof(struct part)) !=
0)
error("Error while allocating memory for SPH particles");
bzero(*parts, *Ngas * sizeof(struct part));
/* Allocate memory to store star particles */
*Nstars = N[STAR];
if (posix_memalign((void*)sparts, spart_align,
*Nstars * sizeof(struct spart)) != 0)
error("Error while allocating memory for star particles");
bzero(*sparts, *Nstars * sizeof(struct spart));
/* Allocate memory to store all particles */
Ndm = N[1];
*Ngparts = N[1] + N[0];
Ndm = N[DM];
*Ngparts = N[GAS] + N[DM] + N[STAR];
if (posix_memalign((void*)gparts, gpart_align,
*Ngparts * sizeof(struct gpart)) != 0)
error("Error while allocating memory for gravity particles");
......@@ -491,6 +502,11 @@ void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
darkmatter_read_particles(*gparts, list, &num_fields);
break;
case STAR:
Nparticles = *Nstars;
star_read_particles(*sparts, list, &num_fields);
break;
default:
message("Particle Type %d not yet supported. Particles ignored", ptype);
}
......@@ -507,9 +523,12 @@ void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
/* Prepare the DM particles */
if (!dry_run) prepare_dm_gparts(*gparts, Ndm);
/* Now duplicate the hydro particle into gparts */
/* Duplicate the hydro particles into gparts */
if (!dry_run) duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
/* Duplicate the star particles into gparts */
if (!dry_run) duplicate_star_gparts(*sparts, *gparts, *Nstars, Ndm + *Ngas);
/* message("Done Reading particles..."); */
/* Clean up */
......@@ -541,18 +560,20 @@ void write_output_single(struct engine* e, const char* baseName,
hid_t h_file = 0, h_grp = 0;
const size_t Ngas = e->s->nr_parts;
const size_t Nstars = e->s->nr_sparts;
const size_t Ntot = e->s->nr_gparts;
int periodic = e->s->periodic;
int numFiles = 1;
struct part* parts = e->s->parts;
struct gpart* gparts = e->s->gparts;
struct gpart* dmparts = NULL;
struct spart* sparts = e->s->sparts;
static int outputCount = 0;
/* Number of unassociated gparts */
const size_t Ndm = Ntot > 0 ? Ntot - Ngas : 0;
const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
long long N_total[NUM_PARTICLE_TYPES] = {Ngas, Ndm, 0};
long long N_total[NUM_PARTICLE_TYPES] = {Ngas, Ndm, 0, 0, Nstars, 0};
/* File name */
char fileName[FILENAME_BUFFER_SIZE];
......@@ -729,6 +750,11 @@ void write_output_single(struct engine* e, const char* baseName,
darkmatter_write_particles(dmparts, list, &num_fields);
break;
case STAR:
N = Nstars;
star_write_particles(sparts, list, &num_fields);
break;
default:
error("Particle Type %d not yet supported. Aborting", ptype);
}
......
......@@ -31,7 +31,8 @@
void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
size_t* Ngas, size_t* Ndm, int* periodic, int* flag_entropy,
struct spart** sparts, size_t* Ngas, size_t* Ndm,
size_t* Nstars, int* periodic, int* flag_entropy,
int dry_run);
void write_output_single(struct engine* e, const char* baseName,
......
......@@ -1770,8 +1770,9 @@ void space_init_gparts(struct space *s) {
*/
void space_init(struct space *s, const struct swift_params *params,
double dim[3], struct part *parts, struct gpart *gparts,
size_t Npart, size_t Ngpart, int periodic, int gravity,
int verbose, int dry_run) {
struct spart *sparts, size_t Npart, size_t Ngpart,
size_t Nspart, int periodic, int gravity, int verbose,
int dry_run) {
/* Clean-up everything */
bzero(s, sizeof(struct space));
......@@ -1789,6 +1790,9 @@ void space_init(struct space *s, const struct swift_params *params,
s->nr_gparts = Ngpart;
s->size_gparts = Ngpart;
s->gparts = gparts;
s->nr_sparts = Nspart;
s->size_sparts = Nspart;
s->sparts = sparts;
s->nr_queues = 1; /* Temporary value until engine construction */
/* Decide on the minimal top-level cell size */
......@@ -1848,6 +1852,11 @@ void space_init(struct space *s, const struct swift_params *params,
gparts[k].x[1] += shift[1];
gparts[k].x[2] += shift[2];
}
for (size_t k = 0; k < Nspart; k++) {
sparts[k].x[0] += shift[0];
sparts[k].x[1] += shift[1];
sparts[k].x[2] += shift[2];
}
}
if (!dry_run) {
......@@ -1879,9 +1888,23 @@ void space_init(struct space *s, const struct swift_params *params,
if (gparts[k].x[j] < 0 || gparts[k].x[j] >= dim[j])
error("Not all g-particles are within the specified domain.");
}
/* Same for the sparts */
if (periodic) {
for (size_t k = 0; k < Nspart; k++)
for (int j = 0; j < 3; j++) {
while (sparts[k].x[j] < 0) sparts[k].x[j] += dim[j];
while (sparts[k].x[j] >= dim[j]) sparts[k].x[j] -= dim[j];
}
} else {
for (size_t k = 0; k < Nspart; k++)
for (int j = 0; j < 3; j++)
if (sparts[k].x[j] < 0 || sparts[k].x[j] >= dim[j])
error("Not all s-particles are within the specified domain.");
}
}
/* Allocate the extra parts array. */
/* Allocate the extra parts array for the gas particles. */
if (Npart > 0) {
if (posix_memalign((void *)&s->xparts, xpart_align,
Npart * sizeof(struct xpart)) != 0)
......@@ -1893,6 +1916,7 @@ void space_init(struct space *s, const struct swift_params *params,
space_init_parts(s);
space_init_xparts(s);
space_init_gparts(s);
// space_init_sparts(s); // MATTHIEU
/* Init the space lock. */
if (lock_init(&s->lock) != 0) error("Failed to create space spin-lock.");
......
......@@ -108,6 +108,9 @@ struct space {
/*! The total number of g-parts in the space. */
size_t nr_gparts, size_gparts;
/*! The total number of g-parts in the space. */
size_t nr_sparts, size_sparts;
/*! The particle data (cells have pointers to this). */
struct part *parts;
......@@ -117,6 +120,9 @@ struct space {
/*! The g-particle data (cells have pointers to this). */
struct gpart *gparts;
/*! The s-particle data (cells have pointers to this). */
struct spart *sparts;
/*! General-purpose lock for this space. */
swift_lock_type lock;
......@@ -152,8 +158,9 @@ int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
double *shift);
void space_init(struct space *s, const struct swift_params *params,
double dim[3], struct part *parts, struct gpart *gparts,
size_t Npart, size_t Ngpart, int periodic, int gravity,
int verbose, int dry_run);
struct spart *sparts, size_t Npart, size_t Ngpart,
size_t Nspart, int periodic, int gravity, int verbose,
int dry_run);
void space_sanitize(struct space *s);
void space_map_cells_pre(struct space *s, int full,
void (*fun)(struct cell *c, void *data), void *data);
......
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
* 2016 Tom Theuns (tom.theuns@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/>.
*
******************************************************************************/
#ifndef SWIFT_DEFAULT_GRAVITY_H
#define SWIFT_DEFAULT_GRAVITY_H
#include <float.h>
#include "minmax.h"
/**
* @brief Computes the gravity time-step of a given particle due to self-gravity
*
* @param gp Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static float
gravity_compute_timestep_self(const struct gpart* const gp) {
const float ac2 = gp->a_grav[0] * gp->a_grav[0] +
gp->a_grav[1] * gp->a_grav[1] +
gp->a_grav[2] * gp->a_grav[2];
const float ac = (ac2 > 0.f) ? sqrtf(ac2) : FLT_MIN;
const float dt = sqrtf(2.f * const_gravity_eta * gp->epsilon / ac);
return dt;
}
/**
* @brief Initialises the g-particles for the first time
*
* This function is called only once just after the ICs have been
* read in to do some conversions.
*
* @param gp The particle to act upon
*/
__attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
struct gpart* gp) {
gp->ti_begin = 0;
gp->ti_end = 0;
gp->epsilon = 0.; // MATTHIEU
}
/**
* @brief Prepares a g-particle for the gravity calculation
*
* Zeroes all the relevant arrays in preparation for the sums taking place in
* the variaous tasks
*
* @param gp The particle to act upon
*/
__attribute__((always_inline)) INLINE static void gravity_init_gpart(
struct gpart* gp) {
/* Zero the acceleration */
gp->a_grav[0] = 0.f;
gp->a_grav[1] = 0.f;
gp->a_grav[2] = 0.f;
}
/**
* @brief Finishes the gravity calculation.
*
* Multiplies the forces and accelerations by the appropiate constants
*
* @param gp The particle to act upon
* @param const_G Newton's constant in internal units
*/
__attribute__((always_inline)) INLINE static void gravity_end_force(
struct gpart* gp, float const_G) {
/* Let's get physical... */
gp->a_grav[0] *= const_G;
gp->a_grav[1] *= const_G;
gp->a_grav[2] *= const_G;
}
/**
* @brief Kick the additional variables
*
* @param gp The particle to act upon
* @param dt The time-step for this kick
* @param half_dt The half time-step for this kick
*/
__attribute__((always_inline)) INLINE static void gravity_kick_extra(
struct gpart* gp, float dt, float half_dt) {}
#endif /* SWIFT_DEFAULT_GRAVITY_H */
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (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/>.
*
******************************************************************************/
#ifndef SWIFT_DEFAULT_STAR_DEBUG_H
#define SWIFT_DEFAULT_STAR_DEBUG_H
__attribute__((always_inline)) INLINE static void star_debug_particle(
const struct spart* p) {
printf(
"x=[%.3e,%.3e,%.3e], "
"v_full=[%.3e,%.3e,%.3e] \n t_begin=%d, t_end=%d\n",
p->x[0], p->x[1], p->x[2], p->v_full[0], p->v_full[1], p->v_full[2],
p->ti_begin, p->ti_end);
}
#endif /* SWIFT_DEFAULT_STAR_DEBUG_H */
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (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/>.
*
******************************************************************************/
#ifndef SWIFT_DEFAULT_STAR_IO_H
#define SWIFT_DEFAULT_STAR_IO_H
#include "io_properties.h"
/**
* @brief Specifies which s-particle fields to read from a dataset
*
* @param sparts The s-particle array.
* @param list The list of i/o properties to read.
* @param num_fields The number of i/o fields to read.
*/
void star_read_particles(struct spart* sparts, struct io_props* list,
int* num_fields) {
/* Say how much we want to read */
*num_fields = 4;
/* List what we want to read */
list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
UNIT_CONV_LENGTH, sparts, x);
list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
UNIT_CONV_SPEED, sparts, v);
list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
sparts, mass);
list[3] = io_make_input_field("ParticleIDs", LONGLONG, 1, COMPULSORY,
UNIT_CONV_NO_UNITS, sparts, id);
}
/**
* @brief Specifies which s-particle fields to write to a dataset
*
* @param sparts The s-particle array.
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
*/
void star_write_particles(struct spart* sparts, struct io_props* list,
int* num_fields) {
/* Say how much we want to read */
*num_fields = 4;
/* List what we want to read */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
sparts, x);
list[1] =
io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, sparts, v);
list[2] =
io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, sparts, mass);
list[3] = io_make_output_field("ParticleIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS,
sparts, id);
}
#endif /* SWIFT_DEFAULT_STAR_IO_H */
/*******************************************************************************
* 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