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

Move the code for searching particles/cells into a separate utility class #234

Merged
merged 17 commits into from
Jan 13, 2025
Merged

Move the code for searching particles/cells into a separate utility class #234

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
d2fc934
initial BoxSearch implementation (unused)
petercamps Jan 11, 2025
d8cdfb5
use new search class in particle snapshot
petercamps Jan 11, 2025
4a968d9
update comments
petercamps Jan 11, 2025
b5b96cc
use box search in cell snapshot
petercamps Jan 11, 2025
cd65318
use new search in tetra grid
petercamps Jan 11, 2025
1e83b65
remove empty destructor
petercamps Jan 11, 2025
e1a6710
fix tetra contains
petercamps Jan 12, 2025
10e3c6f
add nr of entities to log
petercamps Jan 12, 2025
715e5b2
add docs
petercamps Jan 12, 2025
66faf0d
adjust logging
petercamps Jan 12, 2025
fd63df2
improve for small entity numbers
petercamps Jan 12, 2025
7c329d4
use template argument i.o. std::function<>
petercamps Jan 12, 2025
43c21bb
adjust nr of blocks
petercamps Jan 12, 2025
9ca21a1
replace callbacks by iterators
petercamps Jan 13, 2025
d1f4f66
remove duplicates from entities for ray
petercamps Jan 13, 2025
e47cc9d
restore test for zero entities
petercamps Jan 13, 2025
a771a2f
include functional
petercamps Jan 13, 2025
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
273 changes: 44 additions & 229 deletions SKIRT/core/CellSnapshot.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,215 +13,11 @@

////////////////////////////////////////////////////////////////////

namespace
Box CellSnapshot::boxForCell(int m) const
{
// returns a Box object corresponding to the six elements in the given properties array starting at the given index
Box box(const Array& prop, int i)
{
return Box(prop[i], prop[i + 1], prop[i + 2], prop[i + 3], prop[i + 4], prop[i + 5]);
}

// builds a smart grid in the specified spatial direction (box+0=x, box+1=y, box+2=z) and with the specified size,
// and stores it in output parameter "grid", along with the minimum and maximum coordinate enclosing the cells
void makegrid(const vector<Array>& propv, int dir, int gridsize, Array& grid, double& cmin, double& cmax)
{
int n = propv.size();

// find the spatial range of the cells in the specified direction
cmin = +std::numeric_limits<double>::infinity();
cmax = -std::numeric_limits<double>::infinity();
for (const Array& prop : propv)
{
cmin = min(cmin, prop[dir]);
cmax = max(cmax, prop[dir + 3]);
}

// determine the cell distribution by binning at a decent resolution
int nbins = gridsize * 100;
double binwidth = (cmax - cmin) / nbins;
vector<int> bins(nbins);
for (const Array& prop : propv)
{
double center = 0.5 * (prop[dir] + prop[dir + 3]);
bins[static_cast<int>((center - cmin) / binwidth)] += 1;
}

// determine grid separation points based on the cumulative distribution
grid.resize(gridsize + 1);
grid[0] = -std::numeric_limits<double>::infinity();
int percell = n / gridsize; // target number of particles per cell
int cumul = 0; // cumulative number of particles in processed bins
int gridindex = 1; // index of the next grid separation point to be filled
for (int binindex = 0; binindex < nbins; binindex++)
{
cumul += bins[binindex];
if (cumul > percell * gridindex)
{
grid[gridindex] = cmin + (binindex + 1) * binwidth;
gridindex += 1;
if (gridindex >= gridsize) break;
}
}
grid[gridsize] = std::numeric_limits<double>::infinity();
}

// returns the linear index for element (i,j,k) in a p*p*p table
inline int index(int p, int i, int j, int k)
{
return ((i * p) + j) * p + k;
}
}

////////////////////////////////////////////////////////////////////

// This is a helper class for organizing cuboidal cells in a smart grid, so that
// it is easy to retrieve the first cell that overlaps a given point in space.
// The Box object on which this class is based specifies a cuboid guaranteed to
// enclose all cells in the grid.
class CellSnapshot::CellGrid : public Box
{
// data members initialized during construction
const vector<Array>& _propv; // reference to the original list of cells
int _i; // the box index in the properties list of each cell
int _p; // number of grid cells in each spatial direction
Array _xgrid, _ygrid, _zgrid; // the m+1 grid separation points for each spatial direction
vector<vector<int>> _listv; // the m*m*m lists of indices for cells overlapping each grid cell
int _pmin, _pmax, _ptotal; // minimum, maximum nr of cells in list; total nr of cells in listv

public:
// The constructor creates a cuboidal grid of the specified number of grid cells in each
// spatial direction, and for each of the grid cells it builds a list of all cells
// overlapping the grid cell. In an attempt to distribute the cells evenly over the
// grid cells, the sizes of the grid cells in each spatial direction are chosen so that
// the cell centers are evenly distributed over the grid cells.
CellGrid(const vector<Array>& propv, int boxindex, int gridsize) : _propv(propv), _i(boxindex), _p(gridsize)
{
// build the grids in each spatial direction
double xmin, ymin, zmin, xmax, ymax, zmax;
makegrid(propv, boxindex + 0, gridsize, _xgrid, xmin, xmax);
makegrid(propv, boxindex + 1, gridsize, _ygrid, ymin, ymax);
makegrid(propv, boxindex + 2, gridsize, _zgrid, zmin, zmax);
setExtent(xmin, ymin, zmin, xmax, ymax, zmax);

// make room for p*p*p grid cells
_listv.resize(gridsize * gridsize * gridsize);

// add each cell to the list for every grid cell that it overlaps
int n = propv.size();
for (int m = 0; m != n; ++m)
{
Box cell = box(propv[m], boxindex);

// find indices for first and last grid cell overlapped by cell, in each spatial direction
int i1 = NR::locateClip(_xgrid, cell.xmin());
int i2 = NR::locateClip(_xgrid, cell.xmax());
int j1 = NR::locateClip(_ygrid, cell.ymin());
int j2 = NR::locateClip(_ygrid, cell.ymax());
int k1 = NR::locateClip(_zgrid, cell.zmin());
int k2 = NR::locateClip(_zgrid, cell.zmax());

// add the cell to all grid cells in that 3D range
for (int i = i1; i <= i2; i++)
for (int j = j1; j <= j2; j++)
for (int k = k1; k <= k2; k++)
{
_listv[index(gridsize, i, j, k)].push_back(m);
}
}

// calculate statistics
_pmin = n;
_pmax = 0;
_ptotal = 0;
for (int index = 0; index < gridsize * gridsize * gridsize; index++)
{
int size = _listv[index].size();
_pmin = min(_pmin, size);
_pmax = max(_pmax, size);
_ptotal += size;
}
}

// This function returns the smallest number of cells overlapping a single grid cell.
int minCellRefsPerCell() const { return _pmin; }

// This function returns the largest number of cells overlapping a single grid cell.
int maxCellRefsPerCell() const { return _pmax; }

// This function returns the total number of cell references for all cells in the grid.
int totalCellRefs() const { return _ptotal; }

// This function returns the index (in the list originally passed to the constructor)
// of the first cell in the list that overlaps the specified position,
// or -1 if none of the cells in the list overlap the specified position.
int cellIndexFor(Position r) const
{
// locate the grid cell containing the specified position
int i = NR::locateClip(_xgrid, r.x());
int j = NR::locateClip(_ygrid, r.y());
int k = NR::locateClip(_zgrid, r.z());

// search the list of cells for that grid cell
for (int m : _listv[index(_p, i, j, k)])
{
if (box(_propv[m], _i).contains(r)) return m;
}
return -1;
}

// This function replaces the contents of the specified entity collection by the set of cells
// that overlap the path with specified starting point and direction.
// The weight of a cell is given by the length of the path segment inside the cell.
void getEntities(EntityCollection& entities, Position bfr, Direction bfk) const
{
// use the values in these variables only after an intersection call that returns true
double smin, smax;

// initialize the output collection
entities.clear();

// verify that the path intersects the domain
if (extent().intersects(bfr, bfk, smin, smax))
{
// find the indices for first and last grid cell, in each spatial direction,
// overlapped by the bounding box of the path's intersection with the domain
Box pathbox(bfr + smin * bfk, bfr + smax * bfk);
int i1 = NR::locateClip(_xgrid, pathbox.xmin());
int i2 = NR::locateClip(_xgrid, pathbox.xmax());
int j1 = NR::locateClip(_ygrid, pathbox.ymin());
int j2 = NR::locateClip(_ygrid, pathbox.ymax());
int k1 = NR::locateClip(_zgrid, pathbox.zmin());
int k2 = NR::locateClip(_zgrid, pathbox.zmax());
if (i1 > i2) std::swap(i1, i2); // fix reversed coords for negative bfk components
if (j1 > j2) std::swap(j1, j2);
if (k1 > k2) std::swap(k1, k2);

// loop over all grid cells in that 3D range
for (int i = i1; i <= i2; i++)
for (int j = j1; j <= j2; j++)
for (int k = k1; k <= k2; k++)
{
// if the path intersects the grid cell
Box gridcellbox(_xgrid[i], _ygrid[j], _zgrid[k], _xgrid[i + 1], _ygrid[j + 1], _zgrid[k + 1]);
if (gridcellbox.intersects(bfr, bfk, smin, smax))
{
// loop over all cells in this grid cell
for (int m : _listv[index(_p, i, j, k)])
{
// if the path intersects the cell, add the cell to the output collection
if (box(_propv[m], _i).intersects(bfr, bfk, smin, smax)) entities.add(m, smax - smin);
}
}
}
}
}
};

////////////////////////////////////////////////////////////////////

CellSnapshot::~CellSnapshot()
{
delete _grid;
const auto& prop = _propv[m];
int i = boxIndex();
return Box(prop[i], prop[i + 1], prop[i + 2], prop[i + 3], prop[i + 4], prop[i + 5]);
}

////////////////////////////////////////////////////////////////////
Expand Down Expand Up @@ -264,14 +60,14 @@ void CellSnapshot::readAndClose()
if (maxT && prop[temperatureIndex()] > maxT)
numIgnored++;
else
originalMass = max(0., massIndex() >= 0 ? prop[massIndex()]
: prop[densityIndex()] * box(prop, boxIndex()).volume());
originalMass =
max(0., massIndex() >= 0 ? prop[massIndex()] : prop[densityIndex()] * boxForCell(m).volume());

double metallicMass = originalMass * (useMetallicity() ? prop[metallicityIndex()] : 1.);
double effectiveMass = metallicMass * multiplier();

Mv[m] = effectiveMass;
_rhov[m] = effectiveMass / box(prop, boxIndex()).volume();
_rhov[m] = effectiveMass / boxForCell(m).volume();

totalOriginalMass += originalMass;
totalMetallicMass += metallicMass;
Expand All @@ -288,17 +84,20 @@ void CellSnapshot::readAndClose()
if (n) NR::cdf(_cumrhov, Mv);
}

// if needed, construct a 3D-grid over the domain, and create a list of cells that overlap each grid cell
// if needed, construct a search structure for the cells
if (hasMassDensityPolicy() || needGetEntities())
{
int gridsize = max(20, static_cast<int>(pow(_propv.size(), 1. / 3.) / 5));
string size = std::to_string(gridsize);
log()->info("Constructing intermediate " + size + "x" + size + "x" + size + " grid for cells...");
_grid = new CellGrid(_propv, boxIndex(), gridsize);
log()->info(" Smallest number of cells per grid cell: " + std::to_string(_grid->minCellRefsPerCell()));
log()->info(" Largest number of cells per grid cell: " + std::to_string(_grid->maxCellRefsPerCell()));
log()->info(" Average number of cells per grid cell: "
+ StringUtils::toString(_grid->totalCellRefs() / double(gridsize * gridsize * gridsize), 'f', 1));
log()->info("Constructing search grid for " + std::to_string(_propv.size()) + " cells...");
auto bounds = [this](int m) { return boxForCell(m); };
auto intersects = [this](int m, const Box& box) { return box.intersects(boxForCell(m)); };
_search.loadEntities(_propv.size(), bounds, intersects);

int nb = _search.numBlocks();
log()->info(" Number of blocks in grid: " + std::to_string(nb * nb * nb) + " (" + std::to_string(nb) + "^3)");
log()->info(" Smallest number of cells per block: " + std::to_string(_search.minEntitiesPerBlock()));
log()->info(" Largest number of cells per block: " + std::to_string(_search.maxEntitiesPerBlock()));
log()->info(" Average number of cells per block: "
+ StringUtils::toString(_search.avgEntitiesPerBlock(), 'f', 1));
}
}

Expand All @@ -309,8 +108,8 @@ Box CellSnapshot::extent() const
// if there are no cells, return an empty box
if (_propv.empty()) return Box();

// if there is a cell grid, ask it to return the extent (it is already calculated)
if (_grid) return _grid->extent();
// if there is a search structure, ask it to return the extent (it is already calculated)
if (_search.numBlocks()) return _search.extent();

// otherwise find the spatial range of the cells
double xmin = +std::numeric_limits<double>::infinity();
Expand Down Expand Up @@ -342,7 +141,7 @@ int CellSnapshot::numEntities() const

double CellSnapshot::volume(int m) const
{
return box(_propv[m], boxIndex()).volume();
return boxForCell(m).volume();
}

////////////////////////////////////////////////////////////////////
Expand All @@ -356,8 +155,11 @@ double CellSnapshot::density(int m) const

double CellSnapshot::density(Position bfr) const
{
int m = _grid ? _grid->cellIndexFor(bfr) : -1;
return m >= 0 ? _rhov[m] : 0.;
for (int m : _search.entitiesFor(bfr))
{
if (boxForCell(m).contains(bfr)) return _rhov[m];
}
return 0.;
}

////////////////////////////////////////////////////////////////////
Expand All @@ -371,14 +173,14 @@ double CellSnapshot::mass() const

Position CellSnapshot::position(int m) const
{
return Position(box(_propv[m], boxIndex()).center());
return Position(boxForCell(m).center());
}

////////////////////////////////////////////////////////////////////

Position CellSnapshot::generatePosition(int m) const
{
return random()->position(box(_propv[m], boxIndex()));
return random()->position(boxForCell(m));
}

////////////////////////////////////////////////////////////////////
Expand All @@ -405,14 +207,27 @@ const Array& CellSnapshot::properties(int m) const

void CellSnapshot::getEntities(EntityCollection& entities, Position bfr) const
{
entities.addSingle(_grid->cellIndexFor(bfr));
for (int m : _search.entitiesFor(bfr))
{
if (boxForCell(m).contains(bfr))
{
entities.addSingle(m);
return;
}
}
entities.clear();
}

////////////////////////////////////////////////////////////////////

void CellSnapshot::getEntities(EntityCollection& entities, Position bfr, Direction bfk) const
{
_grid->getEntities(entities, bfr, bfk);
entities.clear();
for (int m : _search.entitiesFor(bfr, bfk))
{
double smin, smax;
if (boxForCell(m).intersects(bfr, bfk, smin, smax)) entities.add(m, smax - smin);
};
}

////////////////////////////////////////////////////////////////////
23 changes: 12 additions & 11 deletions SKIRT/core/CellSnapshot.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
#ifndef CELLSNAPSHOT_HPP
#define CELLSNAPSHOT_HPP

#include "BoxSearch.hpp"
#include "Snapshot.hpp"

////////////////////////////////////////////////////////////////////
Expand Down Expand Up @@ -36,10 +37,6 @@ class CellSnapshot : public Snapshot
{
//================= Construction - Destruction =================

public:
/** The destructor deletes the search data structure, if it was constructed. */
~CellSnapshot();

//========== Reading ==========

public:
Expand Down Expand Up @@ -123,20 +120,24 @@ class CellSnapshot : public Snapshot
data structures were not created, invoking this function causes undefined behavior. */
void getEntities(EntityCollection& entities, Position bfr, Direction bfk) const override;

//=================== Private helper functions =====================

private:
/** This function returns the bounding box representing the cell with the given index. If the
index is out of range, the behavior is undefined. */
Box boxForCell(int m) const;

//======================== Data Members ========================

private:
// data members initialized when reading the input file
vector<Array> _propv; // cell properties as imported

// data members initialized after reading the input file if a density policy has been set
Array _rhov; // density for each cell (not normalized)
Array _cumrhov; // normalized cumulative density distribution for cells
double _mass{0.}; // total effective mass

// data members initialized after reading the input file if a density policy has been set
class CellGrid;
CellGrid* _grid{nullptr}; // smart grid for locating the cell at a given location
Array _rhov; // density for each cell (not normalized)
Array _cumrhov; // normalized cumulative density distribution for cells
double _mass{0.}; // total effective mass
BoxSearch _search; // search structure for locating cells
};

////////////////////////////////////////////////////////////////////
Expand Down
Loading