Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

ParticleID #452

Open
wants to merge 17 commits into
base: develop
Choose a base branch
from
Open

ParticleID #452

Show file tree
Hide file tree
Changes from 16 commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
1c2c215
Support domain exchange of 1-byte element arrays
sekelle Mar 9, 2024
8bbb63b
use std::exclusive_scan to avoid overflow
sekelle May 5, 2024
22b4f47
add setMacRadius to updateExpansionCenters
sekelle Apr 11, 2024
eae310a
ParticleData buffer size optimization and harmonization
sekelle May 20, 2024
ccdcd3e
GroupView indexing for SPH and gravity
sekelle May 24, 2024
d7b0e5e
Hierarchical block time-step propagator
sekelle May 24, 2024
b0c29fa
add searchExtFactor to OctreeNsView
sekelle Apr 29, 2024
93a3d7e
make domain peer and halo searches more tolerant for block-time steps
sekelle Apr 29, 2024
3fd7769
warp-coalesced integration functions
sekelle Apr 30, 2024
84f216c
centralize allocGrowthRates
sekelle May 19, 2024
cc4a455
added resizeAcc for bdt prop
sekelle May 20, 2024
ea18537
clang-formatting
sekelle May 25, 2024
866746b
add id field which is conserved to keep track of particles
OsmanSeckinSimsek May 28, 2024
9100ce8
fix compile error in gpu part-data
OsmanSeckinSimsek Jun 12, 2024
fad56a9
Merge remote-tracking branch 'origin/develop' into particleID
OsmanSeckinSimsek Aug 27, 2024
5cd73a7
fix formatting
OsmanSeckinSimsek Aug 27, 2024
de479f4
fix review comments
OsmanSeckinSimsek Oct 8, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion domain/include/cstone/domain/domain.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -264,7 +264,10 @@ class Domain
gatherArrays(sorter.gatherFunc(), sorter.getMap() + global_.numSendDown(), global_.numAssigned(), exchangeStart,
0, std::tie(x, y, z, h, m), util::reverse(scratch));

float invThetaEff = invThetaVecMac(theta_);
// invThetaEff needs to be more strict than the worst case vec mac, because if the tree is reused multiple
// times, an expansion center might have moved outside the cell, which can mark remote cells as halos on ranks
// that won't be found by findPeers
float invThetaEff = invThetaVecMac(theta_) * haloSearchExt_;
std::vector<int> peers = findPeersMac(myRank_, global_.assignment(), global_.octree(), box(), invThetaEff);

if (firstCall_)
Expand Down
2 changes: 2 additions & 0 deletions main/src/init/evrard_init.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -179,6 +179,8 @@ class EvrardGlassSphere : public ISimInitializer<Dataset>

initEvrardFields(d, settings_);

generateParticleIDs(d, rank, numRanks);

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

generateParticleID just initializes an additional field. It can be called in initEvrardFields where all the other fields are initialized as well. Same for the other test cases.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this requires passing rank and numRanks to the init*Fields functions. Should I pass these parameters as well?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're already using the MPI communicator in generateParticleIDs. Might as well use it to obtain rank and numRanks as well. You can take the communicator as an argument to generateParticleIDs and pass MPI_COMM_WORLD at the call site. (One less function to remove MPI_COMM_WORLD from later).

return globalBox;
}

Expand Down
2 changes: 2 additions & 0 deletions main/src/init/gresho_chan.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -157,6 +157,8 @@ class GreshoChan : public ISimInitializer<Dataset>
T massPart = globalBox.lx() * globalBox.ly() * globalBox.lz() * settings_.at("rho") / d.numParticlesGlobal;
initGreshoChanFields(d, settings_, massPart);

generateParticleIDs(d, rank, numRanks);

return globalBox;
}

Expand Down
2 changes: 2 additions & 0 deletions main/src/init/isobaric_cube_init.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -208,6 +208,8 @@ class IsobaricCubeGlass : public ISimInitializer<Dataset>

initIsobaricCubeFields(d, settings_, massPart);

generateParticleIDs(d, rank, numRanks);

return globalBox;
}

Expand Down
2 changes: 2 additions & 0 deletions main/src/init/kelvin_helmholtz_init.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -199,6 +199,8 @@ class KelvinHelmholtzGlass : public ISimInitializer<Dataset>

initKelvinHelmholtzFields(d, settings_, particleMass);

generateParticleIDs(d, rank, numRanks);

return globalBox;
}

Expand Down
2 changes: 2 additions & 0 deletions main/src/init/noh_init.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -145,6 +145,8 @@ class NohGlassSphere : public ISimInitializer<Dataset>

initNohFields(d, settings_);

generateParticleIDs(d, rank, numRanks);

return globalBox;
}

Expand Down
4 changes: 4 additions & 0 deletions main/src/init/sedov_init.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -126,6 +126,8 @@ class SedovGrid : public ISimInitializer<Dataset>

initSedovFields(d, settings_);

generateParticleIDs(d, rank, numRanks);

return globalBox;
}

Expand Down Expand Up @@ -183,6 +185,8 @@ class SedovGlass : public ISimInitializer<Dataset>

initSedovFields(d, settings_);

generateParticleIDs(d, rank, numRanks);

return globalBox;
}

Expand Down
2 changes: 2 additions & 0 deletions main/src/init/turbulence_init.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -139,6 +139,8 @@ class TurbulenceGlass : public ISimInitializer<Dataset>

initTurbulenceHydroFields(d, settings_);

generateParticleIDs(d, rank, numRanks);

return globalBox;
}

Expand Down
23 changes: 23 additions & 0 deletions main/src/init/utils.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -120,6 +120,29 @@ void readFileAttributes(InitSettings& settings, const std::string& settingsFile,
}
}

//! @brief generate particle IDs at the beginning of the simulation initialization
template<class Dataset>
void generateParticleIDs(Dataset& d, int rank, int numRanks)
Copy link
Collaborator

@sekelle sekelle Sep 16, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function doesn't need the full dataset, only a std::span<uint64_t> id . At the call site d.id can be passed, which will either have size 0 if ids are (temporarily) disabled or have size numLocalParticles.

{
std::vector<size_t> ranksLocalParticles(numRanks);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

use uint64_t here to avoid compatibility issues with MacOS

size_t localNumRanks = d.x.size();
// fill ranksLocalParticles with the number of particles per rank
MPI_Allgather(&localNumRanks, 1, MPI_UNSIGNED_LONG, ranksLocalParticles.data(), 1, MPI_UNSIGNED_LONG,
Copy link
Collaborator

@sekelle sekelle Sep 16, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

use MpiType<uint64_t>{} here instead of MPI_UNSIGNED_LONG

MPI_COMM_WORLD);

size_t offset = 0;

for (int i = 0; i < rank; i++)
{
offset += ranksLocalParticles[i];
}

for (size_t i = 0; i < d.x.size(); i++)
{
d.id[i] = offset + i;
}
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These for loops are equivalent to:

std::exclusive_scan(ranksLocalParticles.begin(), ranksLocalParticles.end(), ranksLocalParticles.begin(), uint64_t(0));
std::iota(id.begin(), id.end(), ranksLocalParticles[rank]);

This will work also if id is not actually allocated, which is needed if we want to make it an optional feature later or temporarily disable it by removing id from the list of conserved fields.

}

//! @brief Used to read the default values of dataset attributes
class BuiltinReader
{
Expand Down
2 changes: 2 additions & 0 deletions main/src/init/wind_shock_init.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -220,6 +220,8 @@ class WindShockGlass : public ISimInitializer<Dataset>

initWindShockFields(d, settings_, massPart);

generateParticleIDs(d, rank, numRanks);

return globalBox;
}

Expand Down
2 changes: 1 addition & 1 deletion main/src/propagator/std_hydro.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ class HydroProp : public Propagator<DomainType, DataType>
*
* x, y, z, h and m are automatically considered conserved and must not be specified in this list
*/
using ConservedFields = FieldList<"temp", "vx", "vy", "vz", "x_m1", "y_m1", "z_m1", "du_m1">;
using ConservedFields = FieldList<"temp", "vx", "vy", "vz", "x_m1", "y_m1", "z_m1", "du_m1", "id">;

//! @brief the list of dependent particle fields, these may be used as scratch space during domain sync
using DependentFields =
Expand Down
2 changes: 1 addition & 1 deletion main/src/propagator/std_hydro_grackle.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ class HydroGrackleProp final : public HydroProp<DomainType, DataType>
*
* x, y, z, h and m are automatically considered conserved and must not be specified in this list
*/
using ConservedFields = FieldList<"u", "vx", "vy", "vz", "x_m1", "y_m1", "z_m1", "du_m1">;
using ConservedFields = FieldList<"u", "vx", "vy", "vz", "x_m1", "y_m1", "z_m1", "du_m1", "id">;

//! @brief the list of dependent particle fields, these may be used as scratch space during domain sync
using DependentFields =
Expand Down
2 changes: 1 addition & 1 deletion main/src/propagator/ve_hydro.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ class HydroVeProp : public Propagator<DomainType, DataType>
*
* x, y, z, h and m are automatically considered conserved and must not be specified in this list
*/
using ConservedFields = FieldList<"temp", "vx", "vy", "vz", "x_m1", "y_m1", "z_m1", "du_m1", "alpha">;
using ConservedFields = FieldList<"temp", "vx", "vy", "vz", "x_m1", "y_m1", "z_m1", "du_m1", "alpha", "id">;

//! @brief list of dependent fields, these may be used as scratch space during domain sync
using DependentFields_ =
Expand Down
2 changes: 1 addition & 1 deletion main/src/propagator/ve_hydro_bdt.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,7 @@ class HydroVeBdtProp : public Propagator<DomainType, DataType>
*
* x, y, z, h and m are automatically considered conserved and must not be specified in this list
*/
using ConservedFields = FieldList<"temp", "vx", "vy", "vz", "x_m1", "y_m1", "z_m1", "du_m1", "alpha", "rung">;
using ConservedFields = FieldList<"temp", "vx", "vy", "vz", "x_m1", "y_m1", "z_m1", "du_m1", "alpha", "rung", "id">;

//! @brief list of dependent fields, these may be used as scratch space during domain sync
using DependentFields_ = FieldList<"ax", "ay", "az", "prho", "c", "du", "c11", "c12", "c13", "c22", "c23", "c33",
Expand Down
5 changes: 3 additions & 2 deletions sph/include/sph/particles_data.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -230,6 +230,7 @@ class ParticlesData : public cstone::FieldStates<ParticlesData<AccType>>
FieldVector<unsigned> nc; // number of neighbors of each particle
FieldVector<HydroType> dV11, dV12, dV13, dV22, dV23, dV33; // Velocity gradient components
FieldVector<uint8_t> rung; // rung per particle of previous timestep
FieldVector<uint64_t> id; // unique particle id

//! @brief Indices of neighbors for each particle, length is number of assigned particles * ngmax. CPU version only.
std::vector<cstone::LocalIndex> neighbors;
Expand All @@ -247,7 +248,7 @@ class ParticlesData : public cstone::FieldStates<ParticlesData<AccType>>
"x", "y", "z", "x_m1", "y_m1", "z_m1", "vx", "vy", "vz", "rho", "u", "p",
"prho", "tdpdTrho", "h", "m", "c", "ax", "ay", "az", "du", "du_m1", "c11", "c12",
"c13", "c22", "c23", "c33", "mue", "mui", "temp", "cv", "xm", "kx", "divv", "curlv",
"alpha", "gradh", "keys", "nc", "dV11", "dV12", "dV13", "dV22", "dV23", "dV33", "rung"};
"alpha", "gradh", "keys", "nc", "dV11", "dV12", "dV13", "dV22", "dV23", "dV33", "rung", "id"};

//! @brief dataset prefix to be prepended to fieldNames for structured output
static const inline std::string prefix{};
Expand All @@ -263,7 +264,7 @@ class ParticlesData : public cstone::FieldStates<ParticlesData<AccType>>
{
auto ret = std::tie(x, y, z, x_m1, y_m1, z_m1, vx, vy, vz, rho, u, p, prho, tdpdTrho, h, m, c, ax, ay, az, du,
du_m1, c11, c12, c13, c22, c23, c33, mue, mui, temp, cv, xm, kx, divv, curlv, alpha, gradh,
keys, nc, dV11, dV12, dV13, dV22, dV23, dV33, rung);
keys, nc, dV11, dV12, dV13, dV22, dV23, dV33, rung, id);
#if defined(__clang__) || __GNUC__ > 11
static_assert(std::tuple_size_v<decltype(ret)> == fieldNames.size());
#endif
Expand Down
5 changes: 3 additions & 2 deletions sph/include/sph/particles_data_gpu.cuh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,6 +102,7 @@ public:
DevVector<unsigned> nc; // number of neighbors of each particle
DevVector<HydroType> dV11, dV12, dV13, dV22, dV23, dV33; // Velocity gradient components
DevVector<uint8_t> rung; // rung per particle of previous timestep
DevVector<uint64_t> id; // unique particle id

//! @brief SPH interpolation kernel lookup tables
DevVector<HydroType> wh, whd;
Expand All @@ -118,7 +119,7 @@ public:
"x", "y", "z", "x_m1", "y_m1", "z_m1", "vx", "vy", "vz", "rho", "u", "p",
"prho", "tdpdTrho", "h", "m", "c", "ax", "ay", "az", "du", "du_m1", "c11", "c12",
"c13", "c22", "c23", "c33", "mue", "mui", "temp", "cv", "xm", "kx", "divv", "curlv",
"alpha", "gradh", "keys", "nc", "dV11", "dV12", "dV13", "dV22", "dV23", "dV33", "rung"};
"alpha", "gradh", "keys", "nc", "dV11", "dV12", "dV13", "dV22", "dV23", "dV33", "rung", "id"};

/*! @brief return a tuple of field references
*
Expand All @@ -128,7 +129,7 @@ public:
{
auto ret = std::tie(x, y, z, x_m1, y_m1, z_m1, vx, vy, vz, rho, u, p, prho, tdpdTrho, h, m, c, ax, ay, az, du,
du_m1, c11, c12, c13, c22, c23, c33, mue, mui, temp, cv, xm, kx, divv, curlv, alpha, gradh,
keys, nc, dV11, dV12, dV13, dV22, dV23, dV33, rung);
keys, nc, dV11, dV12, dV13, dV22, dV23, dV33, rung, id);

static_assert(std::tuple_size_v<decltype(ret)> == fieldNames.size());
return ret;
Expand Down
6 changes: 4 additions & 2 deletions sph/include/sph/positions.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,7 @@ void updatePositionsHost(size_t startIndex, size_t endIndex, Dataset& d, const c
template<class Dataset>
void updateTempHost(size_t startIndex, size_t endIndex, Dataset& d)
{
using Tdu = decltype(d.du[0]);
bool haveMui = !d.mui.empty();
auto constCv = idealGasCv(d.muiConst, d.gamma);

Expand All @@ -132,18 +133,19 @@ void updateTempHost(size_t startIndex, size_t endIndex, Dataset& d)
{
auto cv = haveMui ? idealGasCv(d.mui[i], d.gamma) : constCv;
auto u_old = cv * d.temp[i];
d.temp[i] = energyUpdate(u_old, d.minDt, d.minDt_m1, d.du[i], d.du_m1[i]) / cv;
d.temp[i] = energyUpdate(u_old, d.minDt, d.minDt_m1, d.du[i], Tdu(d.du_m1[i])) / cv;
d.du_m1[i] = d.du[i];
}
}

template<class Dataset>
void updateIntEnergyHost(size_t startIndex, size_t endIndex, Dataset& d)
{
using Tdu = decltype(d.du[0]);
#pragma omp parallel for schedule(static)
for (size_t i = startIndex; i < endIndex; i++)
{
d.u[i] = energyUpdate(d.u[i], d.minDt, d.minDt_m1, d.du[i], d.du_m1[i]);
d.u[i] = energyUpdate(d.u[i], d.minDt, d.minDt_m1, d.du[i], Tdu(d.du_m1[i]));
d.du_m1[i] = d.du[i];
}
}
Expand Down
Loading