Commit 35f2cd85 authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge branch 'GravityParticles' into 'master'

Gravity particles

Ok, that's another big one...  This includes all the changes by @tt and @jregan regarding the use of an external potential.

The main changes are:
 - Update to the `unit` with some more const statements.
 - Addition of a set of physical constants and of an object to convert them in the chosen system of units.
 - A new external gravity task and its dependencies.
 - This task calls functions in a new file `potential.h` where users can define their external potentials.
 - A new test case that just makes a bunch of particles orbit a point mass.

See merge request !143
parents bc684e61 35673704
# Define the system of units to use internally.
UnitSystem:
UnitMass_in_cgs: 1.9885e33 # Grams
UnitLength_in_cgs: 3.0856776e21 # Centimeters
UnitVelocity_in_cgs: 1e5 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters for the task scheduling
Scheduler:
nr_queues: 0 # The number of task queues to use. Use 0 to let the system decide.
cell_max_size: 8000000 # Maximal number of interactions per task (this is the default value).
cell_sub_size: 8000000 # Maximal number of interactions per sub-task (this is the default value).
cell_split_size: 400 # Maximal number of particles per cell (this is the default value).
# Parameters governing the time integration
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-3 # The maximal time-step size of the simulation (in internal units).
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2349 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
delta_neighbours: 1. # The tolerance for the targetted number of neighbours.
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
max_ghost_iterations: 30 # Maximal number of iterations allowed to converge towards the smoothing length.
max_smoothing_length: 10. # Maximal smoothing length allowed (in internal units).
# Parameters related to the initial conditions
InitialConditions:
file_name: Sphere.hdf5 # The file to read
h_scaling: 1. # A scaling factor to apply to all smoothing lengths in the ICs.
shift_x: 50. # A shift to apply to all particles read from the ICs (in internal units).
shift_y: 50.
shift_z: 50.
# Parameters govering domain decomposition
DomainDecomposition:
initial_type: m # The initial strategy ("g", "m", "w", or "v"). See documentation for details.
initial_grid_x: 10 # Grid size if the 'g' strategy is chosen.
initial_grid_y: 10
initial_grid_z: 10
repartition_type: b # The re-decomposition strategy ("n", "b", "v", "e" or "x"). See documentation for details.
# External potential parameters
PointMass:
position_x: 50. # location of external point mass in internal units
position_y: 50.
position_z: 50.
mass: 1e10 # mass of external point mass in internal units
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
# 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/>.
#
##############################################################################
import h5py
import sys
import numpy
import math
import random
# Generates a random distriution of particles, for motion in an external potnetial centred at (0,0,0)
# physical constants in cgs
NEWTON_GRAVITY_CGS = 6.672e-8
SOLAR_MASS_IN_CGS = 1.9885e33
PARSEC_IN_CGS = 3.0856776e18
# choice of units
const_unit_length_in_cgs = (1000*PARSEC_IN_CGS)
const_unit_mass_in_cgs = (SOLAR_MASS_IN_CGS)
const_unit_velocity_in_cgs = (1e5)
print "UnitMass_in_cgs: ", const_unit_mass_in_cgs
print "UnitLength_in_cgs: ", const_unit_length_in_cgs
print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs
# derived units
const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
const_G = ((NEWTON_GRAVITY_CGS*const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/(const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs)))
print 'G=', const_G
# Parameters
periodic= 1 # 1 For periodic box
boxSize = 100. #
Radius = boxSize / 4. # maximum radius of particles
G = const_G
Mass = 1e10
N = int(sys.argv[1]) # Number of particles
L = N**(1./3.)
# these are not used but necessary for I/O
rho = 2. # Density
P = 1. # Pressure
gamma = 5./3. # Gas adiabatic index
fileName = "Sphere.hdf5"
#---------------------------------------------------
numPart = N
mass = 1
internalEnergy = P / ((gamma - 1.)*rho)
#--------------------------------------------------
#File
file = h5py.File(fileName, 'w')
# Header
grp = file.create_group("/Header")
grp.attrs["BoxSize"] = boxSize
grp.attrs["NumPart_Total"] = [0, numPart, 0, 0, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [0, numPart, 0, 0, 0, 0]
grp.attrs["Time"] = 0.0
grp.attrs["NumFilesPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
#Particle group
#grp0 = file.create_group("/PartType0")
grp1 = file.create_group("/PartType1")
#generate particle positions
radius = Radius * (numpy.random.rand(N))**(1./3.)
ctheta = -1. + 2 * numpy.random.rand(N)
stheta = numpy.sqrt(1.-ctheta**2)
phi = 2 * math.pi * numpy.random.rand(N)
r = numpy.zeros((numPart, 3))
# r[:,0] = radius * stheta * numpy.cos(phi)
# r[:,1] = radius * stheta * numpy.sin(phi)
# r[:,2] = radius * ctheta
r[:,0] = radius
#
speed = numpy.sqrt(G * Mass / radius)
v = numpy.zeros((numPart, 3))
omega = speed / radius
period = 2.*math.pi/omega
print 'period = minimum = ',min(period), ' maximum = ',max(period)
v[:,0] = -omega * r[:,1]
v[:,1] = omega * r[:,0]
ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
ds[()] = v
v = numpy.zeros(1)
m = numpy.full((numPart, ), mass)
ds = grp1.create_dataset('Masses', (numPart,), 'f')
ds[()] = m
m = numpy.zeros(1)
h = numpy.full((numPart, ), 1.1255 * boxSize / L)
ds = grp1.create_dataset('SmoothingLength', (numPart,), 'f')
ds[()] = h
h = numpy.zeros(1)
u = numpy.full((numPart, ), internalEnergy)
ds = grp1.create_dataset('InternalEnergy', (numPart,), 'f')
ds[()] = u
u = numpy.zeros(1)
ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False)
ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L')
ds[()] = ids
ds = grp1.create_dataset('Coordinates', (numPart, 3), 'd')
ds[()] = r
file.close()
#!/bin/bash
# Generate the initial conditions if they are not present.
if [ ! -e uniformBox.hdf5 ]
then
echo "Generating initial conditions for the point mass potential box example..."
python makeIC.py 10000
fi
../swift -g -t 2 externalPointMass.yml
;
; test energy / angular momentum conservation of test problem
;
@physunits
indir = '/gpfs/data/tt/Codes/Swift-git/swiftsim/examples/'
basefile = 'output_'
nfiles = 657
nfollow = 100 ; number of particles to follow
eout = fltarr(nfollow, nfiles)
ekin = fltarr(nfollow, nfiles)
epot = fltarr(nfollow, nfiles)
tout = fltarr(nfiles)
; set properties of potential
uL = 1e3 * phys.pc ; unit of length
uM = phys.msun ; unit of mass
uV = 1d5 ; unit of velocity
; derived units
constG = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
pcentre = [50.,50.,50.] * 1d3 * pc / uL
mextern = 1d10 * msun / uM
;
;
;
ifile = 0
for ifile=0,nfiles-1 do begin
;for ifile=0,3 do begin
inf = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
time = h5ra(inf, 'Header','Time')
p = h5rd(inf,'PartType1/Coordinates')
v = h5rd(inf,'PartType1/Velocities')
id = h5rd(inf,'PartType1/ParticleIDs')
indx = sort(id)
;
id = id[indx]
for ic=0,2 do begin
tmp = reform(p[ic,*]) & p[ic,*] = tmp[indx]
tmp = reform(v[ic,*]) & v[ic,*] = tmp[indx]
endfor
; calculate energy
dd = size(p,/dimen) & npart = dd[1]
ener = fltarr(npart)
dr = fltarr(npart) & dv = dr
for ic=0,2 do dr[*] = dr[*] + (p[ic,*]-pcentre[ic])^2
for ic=0,2 do dv[*] = dv[*] + v[ic,*]^2
dr = sqrt(dr)
; print,'time = ',time,p[0,0],v[0,0],id[0]
ek = 0.5 * dv
ep = - constG * mextern / dr
ener = ek + ep
tout(ifile) = time
eout(*,ifile) = ener[0:nfollow-1]
ekin(*,ifile) = ek[0:nfollow-1]
epot(*,ifile) = ep[0:nfollow-1]
endfor
; calculate relative energy change
de = 0.0 * eout
for ifile=1, nfiles -1 do de[*,ifile] = (eout[*,ifile]-eout[*,0])/eout[*,0]
end
......@@ -3,6 +3,9 @@
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
* Matthieu Schaller (matthieu.schaller@durham.ac.uk)
* 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
* Angus Lepper (angus.lepper@ed.ac.uk)
* 2016 John A. Regan (john.a.regan@durham.ac.uk)
* 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
......@@ -268,17 +271,6 @@ int main(int argc, char *argv[]) {
MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
#endif
/* Initialize unit system */
struct UnitSystem us;
units_init(&us, params);
if (myrank == 0) {
message("Unit system: U_M = %e g.", us.UnitMass_in_cgs);
message("Unit system: U_L = %e cm.", us.UnitLength_in_cgs);
message("Unit system: U_t = %e s.", us.UnitTime_in_cgs);
message("Unit system: U_I = %e A.", us.UnitCurrent_in_cgs);
message("Unit system: U_T = %e K.", us.UnitTemperature_in_cgs);
}
/* Prepare the domain decomposition scheme */
#ifdef WITH_MPI
struct partition initial_partition;
......@@ -296,6 +288,25 @@ int main(int argc, char *argv[]) {
}
#endif
/* Initialize unit system and constants */
struct UnitSystem us;
struct phys_const prog_const;
units_init(&us, params);
phys_const_init(&us, &prog_const);
if (myrank == 0 && verbose > 0) {
message("Unit system: U_M = %e g.", us.UnitMass_in_cgs);
message("Unit system: U_L = %e cm.", us.UnitLength_in_cgs);
message("Unit system: U_t = %e s.", us.UnitTime_in_cgs);
message("Unit system: U_I = %e A.", us.UnitCurrent_in_cgs);
message("Unit system: U_T = %e K.", us.UnitTemperature_in_cgs);
phys_const_print(&prog_const);
}
/* Initialise the external potential properties */
struct external_potential potential;
if (with_external_gravity) potential_init(params, &us, &potential);
if (with_external_gravity && myrank == 0) potential_print(&potential);
/* Read particles and space information from (GADGET) ICs */
char ICfileName[200] = "";
parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
......@@ -333,6 +344,12 @@ int main(int argc, char *argv[]) {
for (size_t k = 0; k < Ngas; ++k) parts[k].gpart = NULL;
Ngpart = 0;
}
if (!with_hydro) {
free(parts);
parts = NULL;
for (size_t k = 0; k < Ngpart; ++k) if(gparts[k].id > 0) error("Linking problem");
Ngas = 0;
}
/* Get the total number of particles across all nodes. */
long long N_total[2] = {0, 0};
......@@ -396,7 +413,7 @@ int main(int argc, char *argv[]) {
if (myrank == 0) clocks_gettime(&tic);
struct engine e;
engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, engine_policies,
talking);
talking, &prog_const, &potential);
if (myrank == 0) {
clocks_gettime(&toc);
message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
......
......@@ -45,3 +45,11 @@ DomainDecomposition:
initial_grid_z: 10
repartition_type: b # The re-decomposition strategy ("n", "b", "v", "e" or "x"). See documentation for details.
# Parameters related to external potentials
# Point mass external potential
PointMass:
position_x: 50. # location of external point mass in internal units
position_y: 50.
position_z: 50.
mass: 1e10 # mass of external point mass in internal units
......@@ -35,13 +35,15 @@ endif
# List required headers
include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
engine.h swift.h serial_io.h timers.h debug.h scheduler.h proxy.h parallel_io.h \
common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h parser.h
common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h parser.h \
physical_constants.h physical_constants_cgs.h potentials.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
serial_io.c timers.c debug.c scheduler.c proxy.c parallel_io.c \
units.c common_io.c single_io.c multipole.c version.c map.c \
kernel_hydro.c kernel_gravity.c tools.c part.c partition.c clocks.c parser.c
kernel_hydro.c kernel_gravity.c tools.c part.c partition.c clocks.c parser.c \
physical_constants.c potentials.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
......
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
* Matthieu Schaller (matthieu.schaller@durham.ac.uk)
* 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
* 2016 John A. Regan (john.a.regan@durham.ac.uk)
* 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
......
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
* Matthieu Schaller (matthieu.schaller@durham.ac.uk)
* 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
* 2016 John A. Regan (john.a.regan@durham.ac.uk)
* 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
......@@ -122,6 +126,9 @@ struct cell {
/* Tasks for gravity tree. */
struct task *grav_up, *grav_down;
/* Task for external gravity */
struct task *grav_external;
/* Number of tasks that are associated with this cell. */
int nr_tasks;
......
......@@ -70,4 +70,7 @@
#define GADGET2_SPH
//#define DEFAULT_SPH
/* Gravity properties */
#define EXTERNAL_POTENTIAL_POINTMASS
#endif /* SWIFT_CONST_H */
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
* Matthieu Schaller (matthieu.schaller@durham.ac.uk)
* 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
* Angus Lepper (angus.lepper@ed.ac.uk)
* 2016 John A. Regan (john.a.regan@durham.ac.uk)
* 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
......@@ -101,9 +105,12 @@ void engine_make_ghost_tasks(struct engine *e, struct cell *c,
struct cell *super) {
struct scheduler *s = &e->sched;
const int is_with_external_gravity =
(e->policy & engine_policy_external_gravity) ==
engine_policy_external_gravity;
/* Am I the super-cell? */
if (super == NULL && c->nr_tasks > 0) {
if (super == NULL && (c->count > 0 || c->gcount > 0)) {
/* Remember me. */
super = c;
......@@ -111,18 +118,32 @@ void engine_make_ghost_tasks(struct engine *e, struct cell *c,
/* Local tasks only... */
if (c->nodeID == e->nodeID) {
/* Generate the ghost task. */
c->ghost = scheduler_addtask(s, task_type_ghost, task_subtype_none, 0, 0,
c, NULL, 0);
/* Add the drift task. */
c->drift = scheduler_addtask(s, task_type_drift, task_subtype_none, 0, 0,
c, NULL, 0);
/* Add the init task. */
c->init = scheduler_addtask(s, task_type_init, task_subtype_none, 0, 0, c,
NULL, 0);
/* Add the drift task. */
c->drift = scheduler_addtask(s, task_type_drift, task_subtype_none, 0, 0,
c, NULL, 0);
/* Add the kick task. */
c->kick = scheduler_addtask(s, task_type_kick, task_subtype_none, 0, 0, c,
NULL, 0);
if (c->count > 0) {
/* Generate the ghost task. */
c->ghost = scheduler_addtask(s, task_type_ghost, task_subtype_none, 0,
0, c, NULL, 0);
}
if (c->gcount > 0) {
/* Add the external gravity tasks */
if (is_with_external_gravity)
c->grav_external = scheduler_addtask(
s, task_type_grav_external, task_subtype_none, 0, 0, c, NULL, 0);
}
}
}
......@@ -211,27 +232,29 @@ void engine_redistribute(struct engine *e) {
space_parts_sort(s, dest, s->nr_parts, 0, nr_nodes - 1, e->verbose);
/* We need to re-link the gpart partners of parts. */
int current_dest = dest[0];
size_t count_this_dest = 0;
for (size_t k = 0; k < s->nr_parts; ++k) {
if (s->parts[k].gpart != NULL) {
/* As the addresses will be invalidated by the communications, we will */
/* instead store the absolute index from the start of the sub-array */
/* of particles to be sent to a given node. */
/* Recall that gparts without partners have a negative id. */
/* We will restore the pointers on the receiving node later on. */
if (dest[k] != current_dest) {
current_dest = dest[k];
count_this_dest = 0;
}
if (s->nr_parts > 0) {
int current_dest = dest[0];
size_t count_this_dest = 0;
for (size_t k = 0; k < s->nr_parts; ++k) {
if (s->parts[k].gpart != NULL) {
/* As the addresses will be invalidated by the communications, we will */
/* instead store the absolute index from the start of the sub-array */
/* of particles to be sent to a given node. */
/* Recall that gparts without partners have a negative id. */
/* We will restore the pointers on the receiving node later on. */
if (dest[k] != current_dest) {
current_dest = dest[k];
count_this_dest = 0;
}
/* Debug */
/* if(s->parts[k].gpart->id < 0) */
/* error("Trying to link a partnerless gpart !"); */
/* Debug */
/* if(s->parts[k].gpart->id < 0) */
/* error("Trying to link a partnerless gpart !"); */
s->parts[k].gpart->id = count_this_dest;
count_this_dest++;
s->parts[k].gpart->id = count_this_dest;
count_this_dest++;
}
}
}
......@@ -819,7 +842,10 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
ticks tic = getticks();
/* Re-set the proxies. */
for (int k = 0; k < e->nr_proxies; k++) e->proxies[k].nr_parts_out = 0;
for (int k = 0; k < e->nr_proxies; k++) {
e->proxies[k].nr_parts_out = 0;
e->proxies[k].nr_gparts_out = 0;
}
/* Put the parts and gparts into the corresponding proxies. */
for (size_t k = 0; k < *Npart; k++) {
......@@ -839,7 +865,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
/* Re-link the associated gpart with the buffer offset of the part. */
if (s->parts[offset_parts + k].gpart != NULL) {
s->parts[offset_parts + k].gpart->id = e->proxies[pid].nr_parts_in;
s->parts[offset_parts + k].gpart->id = e->proxies[pid].nr_parts_out;
}
/* Load the part and xpart into the proxy. */
......@@ -898,6 +924,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
count_parts_in, count_gparts_in);
}
if (offset_parts + count_parts_in > s->size_parts) {
message("re-allocating parts array.");
s->size_parts = (offset_parts + count_parts_in) * engine_parts_size_grow;
struct part *parts_new = NULL;
struct xpart *xparts_new = NULL;
......@@ -912,8 +939,14 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
free(s->xparts);
s->parts = parts_new;
s->xparts = xparts_new;
for (size_t k = 0; k < offset_parts; k++) {
if (s->parts[k].gpart != NULL) {
s->parts[k].gpart->part = &s->parts[k];
}
}
}
if (offset_gparts + count_gparts_in > s->size_gparts) {
message("re-allocating gparts array.");
s->size_gparts = (offset_gparts + count_gparts_in) * engine_parts_size_grow;
struct gpart *gparts_new = NULL;
if (posix_memalign((void **)&gparts_new, gpart_align,
......@@ -922,6 +955,11 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
memcpy(gparts_new, s->gparts, sizeof(struct gpart) * offset_gparts);
free(s->gparts);
s->gparts = gparts_new;
for (size_t k = 0; k < offset_gparts; k++) {
if (s->gparts[k].id > 0) {
s->gparts[k].part->gpart = &s->gparts[k];
}
}
}
/* Collect the requests for the particle data from the proxies. */
......@@ -976,7 +1014,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
reqs_in[pid + 1] == MPI_REQUEST_NULL &&
reqs_in[pid + 2] == MPI_REQUEST_NULL) {
/* Copy the particle data to the part/xpart/gpart arrays. */
struct proxy *p = &e->proxies[pid >> 1];
struct proxy *p = &e->proxies[pid / 3];
memcpy(&s->parts[offset_parts + count_parts], p->parts_in,
sizeof(struct part) * p->nr_parts_in);
memcpy(&s->xparts[offset_parts + count_parts], p->xparts_in,
......@@ -993,7 +1031,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
for (int k = 0; k < p->nr_gparts_in; k++) {
struct gpart *gp = &s->gparts[offset_gparts + count_gparts + k];
if (gp->id >= 0) {
struct part *p = &s->parts[offset_gparts + count_parts + gp->id];
struct part *p = &s->parts[offset_parts + count_parts + gp->id];
gp->part = p;
p->gpart = gp;
}
......@@ -1004,7 +1042,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
count_gparts += p->nr_gparts_in;