pore_kernels.hpp File Reference

peclet.pnm — the ONE set of stage kernels shared by the single-rank (pore_extraction.hpp) and the distributed (pore_extraction_mpi.hpp) pipelines. More...

#include <algorithm>
#include <array>
#include <cstddef>
#include <cstdint>
#include <Kokkos_Core.hpp>
#include <utility>
#include <vector>

Go to the source code of this file.

Classes
struct	pnm::Pore
struct	pnm::I3
struct	pnm::GridGeo
	Single-rank geometry policy: the field IS the whole periodic grid, owned == field coords, gid == index. More...
struct	pnm::BlockGeo
	Extended-block geometry policy (the MPI pipeline): ghost-inclusive local box <-> global periodic grid. More...
struct	pnm::kernels::FlatFaces< V >
	Per-face field accessors for the film-attachment sweep: label / pair key of face (cell e, d). More...
struct	pnm::kernels::SplitFaces< V >
struct	pnm::kernels::ThroatSlots
	Global unique (patch root, pair key) list -> the throat list ordered by (pair, root fid), the sorted root keys and the (key position -> throat slot) map for the device lookup. More...

Namespaces
namespace	pnm
namespace	pnm::kernels

Typedefs
using	pnm::Exec = Kokkos::DefaultExecutionSpace
using	pnm::Mem = Exec::memory_space
using	pnm::Index = std::int64_t
	Global voxel / face id (== peclet::core::Index; asserted where core is included).
using	pnm::kernels::MD3 = Kokkos::MDRangePolicy< Exec, Kokkos::Rank< 3 > >
using	pnm::kernels::R1 = Kokkos::RangePolicy< Exec >

Functions
template<class T >
void	pnm::uploadVec (const std::vector< T > &h, const Kokkos::View< T *, Mem > &d)
	Bulk host->device upload of a whole std::vector via one deep_copy over an unmanaged host view.
template<class V >
std::vector< typename V::value_type >	pnm::downloadN (const V &d, std::size_t count)
	Download the first count elements of a device view into a fresh std::vector via one deep_copy — only moves what is used.
template<class V >
V::value_type	pnm::readScalar (const V &d)
	Read a device scalar (the changed / counter flags) back to the host.
int	pnm::get_idx (int x, int y, int z, I3 res)
int	pnm::wrapc (int a, int n)
template<class Geo >
MD3	pnm::kernels::ownedRange (const Exec &space, const Geo &geo)
	The owned box of a geometry as an MDRange policy.
template<class Geo , class SdfView >
bool	pnm::kernels::isPeak (const Geo &g, const SdfView &sdf, int lx, int ly, int lz, float cv, Index cgid)
	Local SDF maximum under the global (sdf, gid) lexicographic order (26-neighbourhood).
template<class Geo , class SdfView >
Pore	pnm::kernels::poreAt (const Geo &g, const SdfView &sdf, int lx, int ly, int lz, float ox, float oy, float oz, float sx, float sy, float sz)
	The pore record at a peak: sub-voxel position by the SDF^2-weighted centroid of the 3^3 neighbourhood (positive part only), in global (wrapped) integer coords, radius = sdf there.
template<class Geo , class SdfView >
Index	pnm::kernels::steepest (const Geo &g, const SdfView &sdf, int lx, int ly, int lz, Index cgid, bool ascent)
	Steepest 26-neighbour under the (sdf, gid) lexicographic order: ascent (pores) or descent (solids); returns the neighbour's gid, or cgid itself at a peak.
template<class Geo , class FieldView >
double	pnm::kernels::trilinear (const Geo &g, const FieldView &p, const double gp[3])
	Trilinear interpolation of a cell-centered periodic field at grid position gp (in cells).
bool	pnm::kernels::facesAdjacent (int dx, int dy, int dz, int d, int d2)
	Face adjacency of the per-patch throat CCL: two faces (cell offset (dx,dy,dz), directions d and d2) connect when their doubled center offset satisfies |2*dc + e_d2 - e_d|^2 <= 8 — a generous edge/corner adjacency in pure integer arithmetic (periodic-safe, decomposition- independent).
template<class KeyView >
Index	pnm::kernels::bsearchKey (const KeyView &keys, Index nkeys, typename KeyView::non_const_value_type key)
	Device binary search over a sorted key view; returns the position or -1.
template<class Geo , class SdfView , class PoreView , class CountView >
void	pnm::kernels::detectPores (const Exec &space, const Geo &geo, const SdfView &sdf, std::array< float, 3 > origin, std::array< float, 3 > spacing, const PoreView &pores, const CountView &counter, int maxPores)
	Emit the owned peaks (local SDF maxima, sdf > 0) as Pore records into pores (atomic slot counter counter, capacity maxPores).
template<class Geo , class SdfView , class ParentView >
void	pnm::kernels::initMarkers (const Exec &space, const Geo &geo, const SdfView &sdf, float thr, const ParentView &parent)
	Deep-solid markers: parent(o) = o where sdf < thr, else -1 (owned-dense).
template<class ParentView >
void	pnm::kernels::cclFlatten (const Exec &space, const ParentView &parent, std::size_t n)
	Path compression: every labelled entry points at its root.
template<class Geo , class ParentView , class FlagView >
void	pnm::kernels::cclMergeMarkers (const Exec &space, const Geo &geo, const ParentView &parent, const FlagView &changed)
	One union sweep of the marker CCL (26-connectivity, 13 forward neighbours, union-by-min via atomic_min on the larger root).
template<class ParentView , class FlagView , class MergeFn >
void	pnm::kernels::cclFixpoint (const Exec &space, const ParentView &parent, std::size_t n, const FlagView &changed, const MergeFn &merge)
	Union-find to fixpoint: merge() launches one union sweep (setting changed), then a flatten; repeat until a sweep changes nothing.
template<class Geo , class SdfView , class LabelView , class FlagView >
void	pnm::kernels::floodSweep (const Exec &space, const Geo &geo, const SdfView &sdf, const LabelView &labels, const LabelView &labelsN, const FlagView &changed)
	One Jacobi flood sweep: every unlabelled solid voxel takes the smallest labelled 26-neighbour (reads labels, writes labelsN; the caller swaps).
template<class Geo , class SdfView , class LabelView , class FlagView , class SyncFn >
void	pnm::kernels::floodFixpoint (const Exec &space, const Geo &geo, const SdfView &sdf, LabelView &labels, LabelView &labelsN, const FlagView &changed, const SyncFn &sync)
	Jacobi flood to fixpoint.
template<class SdfView , class EmitFn >
void	pnm::kernels::gradientWalk (const Exec &space, const GridGeo &geo, const SdfView &sdf, bool ascentAll, const EmitFn &emit)
	Single-rank gradient walk: from every owned cell follow the steepest neighbour (ascent for pores, descent for solids; ascentAll = ascent from every cell, the flow-basin variant) for at most 512 steps; emit(o, walker) receives the end point (a gid == field index here).
template<class Geo , class SdfView , class TargetView , class ChaseFn , class FinalFn , class PeakFn >
void	pnm::kernels::forestInit (const Exec &space, const Geo &geo, const SdfView &sdf, const TargetView &target, const ChaseFn &chase, const FinalFn &finalOf, const PeakFn &peakOf)
	Steepest-neighbour forest for the hold-at-ghost root resolution.
template<class Geo , class TargetView >
int	pnm::kernels::resolveHoldAtGhost (const Exec &space, const Geo &geo, const TargetView &target)
	One round of pointer jumping with the hold-at-ghost invariant: chase within the field, adopt a finalized (negative) value, or HOLD at the last in-field cell when the next hop leaves the field (never store a gid the field cannot reach — that strands the chase outside the ghost ring).
template<class Geo , class SegView , class PairView , class CountView >
void	pnm::kernels::boundaryPairs (const Exec &space, const Geo &geo, const SegView &seg, const PairView &pairs, const CountView &cnt, Index maxPairs)
	Boundary pairs: for every owned cell, the (min,max) label pairs across its +x/+y/+z faces where the label changes, appended to the interleaved pairs (atomic cnt, capacity maxPairs).
template<class PairView >
std::vector< std::pair< int, int > >	pnm::kernels::uniquePairs (const PairView &pairs, std::size_t count)
	Interleaved device pairs -> host (l1,l2) list, sorted and unique.
template<class Geo , class SegView >
int	pnm::kernels::maxLabel (const Exec &space, const Geo &geo, const SegView &seg)
	The largest label over the owned box (= the number of pore labels).
template<class Geo , class LabView , class SdfView , class OpenView , class CoreView , class ParentView , class PairView >
void	pnm::kernels::faceInit (const Exec &space, const Geo &geo, const LabView &flowLab, const SdfView &sdf, const OpenView &ox, const OpenView &oy, const OpenView &oz, bool hasOpen, const CoreView &core, const ParentView &parent, const PairView &fpair)
	Interface predicate + CORE tier per owned face: an interface face separates two different positive flow-basin labels with positive openness; CORE = both cells fluid-centered (sdf > 0), which defines the patches by CCL; the rest are FILM faces that attach later.
template<class Geo , class ParentView , class PairView , class FlagView >
void	pnm::kernels::faceMerge (const Exec &space, const Geo &geo, const ParentView &parent, const PairView &fpair, const FlagView &changed)
	One union sweep of the face CCL over the faces with parent >= 0 (same pair key + adjacency; owned-owned adjacency only, see cclMergeMarkers).
template<class Geo , class ParentView , class MinView >
void	pnm::kernels::faceMinGid (const Exec &space, const Geo &geo, const ParentView &parent, const MinView &mg)
	Per-component minimum GLOBAL face id (3*gid+d) of the local face CCL, keyed by the component root: the decomposition-independent patch label.
template<class ParentView , class PairView , class LabView , class RootFn >
void	pnm::kernels::faceLabelCore (const Exec &space, std::size_t nf, const ParentView &parent, const PairView &fpair, const LabView &faceLab, const RootFn &rootOf)
	Core-tier labels: faceLab = rootOf(parent) for the CCL'd faces, kSent for the remaining interface (film) faces, -1 for non-interface faces.
template<class ParentView , class LabView , class RootFn >
void	pnm::kernels::faceLabelLeftover (const Exec &space, std::size_t nf, const ParentView &parent, const LabView &faceLab, const RootFn &rootOf)
	Leftover-tier labels: the still-kSent films (their own CCL tier) get rootOf(parent).
template<class ParentView , class LabView >
void	pnm::kernels::leftoverInit (const Exec &space, std::size_t nf, const LabView &faceLab, const ParentView &parent)
	Leftover tier init: parent = f for the films no core patch reached, -1 otherwise.
template<class Geo , class LabView , class PairView , class CoreView , class NbLab , class NbPair , class FlagView >
void	pnm::kernels::filmAttachSweep (const Exec &space, const Geo &geo, const LabView &faceLab, const PairView &fpair, const CoreView &core, const NbLab &nbLab, const NbPair &nbPair, const FlagView &changed)
	One Jacobi min-propagation sweep of the film attachment: every film face (interface, not core) takes the minimum attached label among its adjacent same-pair faces (nbLab(e, d) / nbPair(e, d) read the neighbour faces — flat single-rank, per-direction exchanged fields on a block).
ThroatSlots	pnm::kernels::throatSlots (std::vector< std::pair< Index, Index > > rootPair)
template<class Geo , class LabView , class FieldView , class OpenView , class FaceLabView , class KeyView , class SlotView , class AncView , class AccView >
void	pnm::kernels::throatFlux (const Exec &space, const Geo &geo, const LabView &flowLab, const FieldView &u, const FieldView &v, const FieldView &w, const OpenView &ox, const OpenView &oy, const OpenView &oz, bool hasOpen, const FaceLabView &faceLab, const KeyView &keyD, const SlotView &slotD, std::size_t nt, const AncView &ancX, const AncView &ancY, const AncView &ancZ, std::array< float, 3 > origin, std::array< float, 3 > spacing, std::array< int, 3 > gdims, const AccView &Q, const AccView &A, const AccView &Cx, const AccView &Cy, const AccView &Cz, const AccView &resid)
	Accumulate over every owned +face where the flow-basin label changes: the openness-weighted MAC flux o·u·A into the per-pore residual (signed) and, for interface faces of a throat patch, into the throat's flow (positive lower -> higher label), open area and area-weighted centroid (min-imaged relative to the LOWER pore's peak-voxel anchor anc{X,Y,Z} — integer-derived, so the periodic-image branch is deterministic under FMA wobble).
std::vector< double >	pnm::kernels::throatPressureDrops (const std::vector< std::pair< int, int > > &throats, const std::vector< Pore > &pores, const std::vector< double > &porePressure, const std::vector< double > &throatArea, const std::vector< double > &Cx, const std::vector< double > &Cy, const std::vector< double > &Cz, const std::vector< double > &ancX, const std::vector< double > &ancY, const std::vector< double > &ancZ, std::array< float, 3 > spacing, std::array< int, 3 > gdims, std::array< double, 3 > grad_p)
	Total-pressure drop per throat: periodic parts + the macroscopic gradient along the throat- anchored two-leg min-image path i -> throat centroid -> j.

Variables
constexpr Index	pnm::kernels::kSent = 0x7ffffffffffffffeLL
	Film-face sentinel of the per-patch throat labelling: "interface face awaiting attachment".

Detailed Description

peclet.pnm — the ONE set of stage kernels shared by the single-rank (pore_extraction.hpp) and the distributed (pore_extraction_mpi.hpp) pipelines.

Every stage kernel is templated on a small GEOMETRY POLICY that maps the owned box a rank iterates over onto the field it reads (the whole periodic grid for a single rank, a ghost-extended block for an MPI rank) and onto GLOBAL voxel ids (the tie-break identity that makes every fixpoint decomposition-independent):

nx, ny, nz owned box (the cells this rank iterates over and writes) g offset of owned (0,0,0) in field coordinates (0 single-rank, the ghost width for a block) cell(lx,ly,lz) field index of FIELD coords (periodic wrap single-rank; plain for a block — ghosts hold the wrapped neighbours) gid(lx,ly,lz) global voxel id of field coords (== cell single-rank) gcoord(lx,ly,lz,...) wrapped global integer coords (the float positions come from these) inOwned(ix,iy,iz) owned coords inside the owned box? (always single-rank: a periodic neighbour IS owned; a block routes cross-block adjacency through its boundary graph instead) owned(ix,iy,iz) dense owned-box index (== cell single-rank) ownedCoords(o,...) the inverse of owned() localOf(gid) field index of a global id or -1 when it has no image in the field cellAtGlobal(gx,gy,gz) field index of (wrapped) global coords

GridGeo is the single-rank policy (trivial: everything is the grid itself); BlockGeo the extended-block policy of the MPI pipeline. The MPI file adds only the halo exchanges, the ownership/merge orchestration and the distributed reductions around these kernels.

Arithmetic is IDENTICAL between the two instantiations by construction (same expression, same evaluation order) — that is what keeps the multi-rank results bit-exact to single-rank.

Definition in file pore_kernels.hpp.