Skip to content

Version for only Gas and BH particles #3

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

Open
robjmcgibbon wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Version for only Gas and BH particles #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
10 changes: 10 additions & 0 deletions compile_reduce_snapshots_Gas_BH.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
#! /bin/bash

# Modules for compiling on cosma
# Using an old version of HDF5 because of the DMantissa9 lossy filter
module purge
module load gnu_comp/14.1.0 openmpi/5.0.3 parallel_hdf5/1.12.3

mpic++ -O3 -Wall -Werror -fopenmp -o reduce_snapshots \
reduce_snapshots_Gas_BH.cpp \
-lhdf5 -lstdc++fs 2>&1 | tee compile_reduce_snapshots_Gas_BH.log
65 changes: 60 additions & 5 deletions reduce_snapshots.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -600,8 +600,8 @@ class SOTable {
std::vector<double> _ZSO;
/*! @brief SO radii. */
std::vector<double> _RSO;
/*! @brief SO masses */
std::vector<double> _MSO;
// /*! @brief SO masses */
// std::vector<double> _MSO; // SO Masses not needed in my version
/*! @brief IDs of the halos. */
std::vector<uint64_t> _haloIDs;

Expand Down Expand Up @@ -1070,6 +1070,7 @@ int main(int argc, char **argv) {
// actual type indices.
const uint_fast32_t typelist[5] = {0, 1, 4, 5, 6};


/// step 2 (from above)

// first partfile read: determine how many files we have (and read some
Expand Down Expand Up @@ -1287,6 +1288,7 @@ int main(int argc, char **argv) {
size_t outcount = 0;
size_t dupcount = 0;
size_t filtercount = 0;
size_t overwrite_count = 0;
// allocate the cell_SOs vector once and then reuse it
// assuming the size does not change much between iterations, this saves
// on reallocations
Expand Down Expand Up @@ -1498,8 +1500,9 @@ int main(int argc, char **argv) {
size_t this_outcount = 0;
size_t this_dupcount = 0;
size_t this_filtercount = 0;
size_t this_overwrite_count =0;
#pragma omp parallel for default(shared) \
reduction(+ : this_keepcount, this_outcount, this_dupcount, this_filtercount)
reduction(+ : this_keepcount, this_outcount, this_dupcount, this_filtercount,this_overwrite_count)
for (int64_t ipart = 0; ipart < size; ++ipart) {
// since we only read the chunks of the particle arrays that are
// actually relevant, the offset stored in the cell meta-data does
Expand Down Expand Up @@ -1611,9 +1614,28 @@ int main(int argc, char **argv) {
// Assign particle to closest halo as function of r/R200c
//
++this_dupcount;
// Define id for old and new halo
const int64_t oldhalo = haloIDs[itype][offset + ipart];
const int64_t newhalo = SOID;
if (d/SOtable.RSO(newhalo) < d/SOtable.RSO(oldhalo)) {
// Get distance from particle to old halo
double dx_old = partpos[3 * totipart] - SOtable.XSO(oldhalo);
if (dx_old < -0.5 * boxSize[0])
dx_old += boxSize[0];
if (dx_old >= 0.5 * boxSize[0])
dx_old -= boxSize[0];
double dy_old = partpos[3 * totipart + 1] - SOtable.YSO(oldhalo);
if (dy_old < -0.5 * boxSize[1])
dy_old += boxSize[1];
if (dy_old >= 0.5 * boxSize[1])
dy_old -= boxSize[1];
double dz_old = partpos[3 * totipart + 2] - SOtable.ZSO(oldhalo);
if (dz_old < -0.5 * boxSize[2])
dz_old += boxSize[2];
if (dz_old >= 0.5 * boxSize[2])
dz_old -= boxSize[2];
double d_old = std::sqrt(dx_old * dx_old + dy_old * dy_old + dz_old * dz_old);
if (d/SOtable.RSO(newhalo) < d_old/SOtable.RSO(oldhalo)) {
++this_overwrite_count;
haloIDs[itype][offset + ipart] = newhalo;
}
}
Expand All @@ -1628,6 +1650,7 @@ int main(int argc, char **argv) {
outcount += this_outcount;
dupcount += this_dupcount;
filtercount += this_filtercount;
overwrite_count += this_overwrite_count;
my_assert(new_sizes[itype][icell] == 0, "Overwriting cell size!");
// register the new number of particles of this type within the top-
// level cell
Expand Down Expand Up @@ -1655,6 +1678,8 @@ int main(int argc, char **argv) {
orphans.size());
timelog(LOGLEVEL_CHUNKLOOPS,
" number of particles filtered out by radius: %zu", filtercount);
timelog(LOGLEVEL_CHUNKLOOPS,
" number of particles with duplicated SO overwritten by r/R200: %zu", overwrite_count);

memory_ifile.add_entry(ifilename.str() + " Orphans", orphans);

Expand Down Expand Up @@ -1700,7 +1725,37 @@ int main(int argc, char **argv) {
#pragma omp critical
++new_sizes[orphan.type][orphan.cell];
++this_keepcount;
}
} else {
// duplicate! particle has been assigned to multiple halo ids
// We must decide which halo to assign this particle too.
//
// Assign particle to closest halo as function of r/R200c
//
//++this_dupcount;
// Define halo id for old and new haloes
const int64_t oldhalo = haloIDs[orphan.type][orphan.index];
const int64_t newhalo = iSO;
// Get distance from particle to old halo
double dx_old = orphan.position[0] - SOtable.XSO(oldhalo);
if (dx_old < -0.5 * boxSize[0])
dx_old += boxSize[0];
if (dx_old >= 0.5 * boxSize[0])
dx_old -= boxSize[0];
double dy_old = orphan.position[1] - SOtable.YSO(oldhalo);
if (dy_old < -0.5 * boxSize[1])
dy_old += boxSize[1];
if (dy_old >= 0.5 * boxSize[1])
dy_old -= boxSize[1];
double dz_old = orphan.position[2] - SOtable.ZSO(oldhalo);
if (dz_old < -0.5 * boxSize[2])
dz_old += boxSize[2];
if (dz_old >= 0.5 * boxSize[2])
dz_old -= boxSize[2];
double d_old = std::sqrt(dx_old * dx_old + dy_old * dy_old + dz_old * dz_old);
if (d/SOtable.RSO(newhalo) < d_old/SOtable.RSO(oldhalo)) {
//++this_overwrite_count;
haloIDs[orphan.type][orphan.index] = newhalo;
}
}
}
}
Expand Down
Loading