diff --git a/CMakeLists.txt b/CMakeLists.txt index 7d792c9..96ae58d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -5,7 +5,7 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) project(plsm - VERSION 1.1.0 + VERSION 2.0.0 LANGUAGES CXX ) @@ -26,7 +26,6 @@ set(PLSM_HEADERS ${PLSM_HEADER_DIR}/refine/PolylineDetector.h ${PLSM_HEADER_DIR}/refine/RegionDetector.h ${PLSM_HEADER_DIR}/CompactFlat.h - ${PLSM_HEADER_DIR}/ContextUtility.h ${PLSM_HEADER_DIR}/EnumIndexed.h ${PLSM_HEADER_DIR}/Interval.h ${PLSM_HEADER_DIR}/IntervalRange.h diff --git a/include/plsm/ContextUtility.h b/include/plsm/ContextUtility.h deleted file mode 100644 index e3da608..0000000 --- a/include/plsm/ContextUtility.h +++ /dev/null @@ -1,103 +0,0 @@ -#pragma once - -#include - -namespace plsm -{ -/*! - * @brief Tag type for host memory space - */ -struct OnHost -{ -}; - -/*! - * @brief Tag type for device memory space - */ -struct OnDevice -{ -}; - -//! Tag instance for host memory space -inline constexpr OnHost onHost{}; -//! Tag instance for device memory space -inline constexpr OnDevice onDevice{}; - -namespace detail -{ -/*! - * @brief Determine the view type appropriate for the given memory space - */ -template -struct ContextualViewTypeHelper -{ -}; - -/*! @cond */ -template -struct ContextualViewTypeHelper -{ - using Type = typename TDualView::t_host; -}; - -template -struct ContextualViewTypeHelper -{ - using Type = typename TDualView::t_dev; -}; -/*! @endcond */ - -/*! - * @brief Determine the part of a Kokkos::DualView appropriate for the given - * memory space - */ -template -using ContextualViewType = - typename ContextualViewTypeHelper::Type; - -//@{ -/*! - * @brief Get the part of a Kokkos::DualView appropriate for the given memory - * space - */ -template -KOKKOS_INLINE_FUNCTION -const ContextualViewType& -getContextualView(const TDualView& dualView, ::plsm::OnHost) -{ - return dualView.h_view; -} - -template -KOKKOS_INLINE_FUNCTION -const ContextualViewType& -getContextualView(const TDualView& dualView, ::plsm::OnDevice) -{ - return dualView.d_view; -} -//@} - -//@{ -/*! - * @brief Syncronize a dual view in the specified direction - */ -template -void -syncUpdate(TDualView& dualView, ::plsm::OnHost) -{ - dualView.modify_device(); - Kokkos::resize(dualView.h_view, dualView.d_view.extent(0)); - dualView.sync_host(); -} - -template -void -syncUpdate(TDualView& dualView, ::plsm::OnDevice) -{ - dualView.modify_host(); - Kokkos::resize(dualView.d_view, dualView.h_view.extent(0)); - dualView.sync_device(); -} -//@} -} // namespace detail -} // namespace plsm diff --git a/include/plsm/Subpaving.h b/include/plsm/Subpaving.h index fdc48c7..eeb2155 100644 --- a/include/plsm/Subpaving.h +++ b/include/plsm/Subpaving.h @@ -4,16 +4,20 @@ #include #include -#include -#include #include +#include #include #include #include namespace plsm { +namespace test +{ +template +struct SubpavingTester; +} /*! * @brief A set of non-overlapping "tiles" which cover the space within a given * Region of an N-dimensional lattice @@ -33,13 +37,20 @@ namespace plsm * @test benchmark_Subpaving.cpp */ template + typename TItemData = IdType, typename TMemSpace = DefaultMemSpace> class Subpaving { + template + friend class ::plsm::test::SubpavingTester; + template friend class detail::Refiner; + static_assert(Kokkos::is_memory_space{}); + public: + //! The memory space for the subpaving data + using MemorySpace = TMemSpace; //! Underlying type for scalar representation for the lattice using ScalarType = TScalar; //! The type to represent lattice regions @@ -53,21 +64,30 @@ class Subpaving //! The subdivision Zone using ZoneType = Zone; - //! The type for the set of zones (host/device synchronized) - using ZonesDualView = Kokkos::DualView; //! The type for the set of zones on the given memory space - template - using ZonesView = detail::ContextualViewType; + using ZonesView = Kokkos::View; + //! Read-only random-access view of zones + using ZonesRAView = + Kokkos::View; //! The subpaving Tile using TileType = Tile; - //! The type for the set of tiles (host/device synchronized) - using TilesDualView = Kokkos::DualView; //! The type for the set of tiles on the given memory space - template - using TilesView = detail::ContextualViewType; + using TilesView = Kokkos::View; + //! Read-only random-access view of tiles + using TilesRAView = + Kokkos::View; - Subpaving() = delete; + using HostMirrorSpace = typename TilesView::traits::host_mirror_space; + using HostMirror = + Subpaving; + +private: + template + friend class ::plsm::Subpaving; + +public: + Subpaving() = default; /*! * @brief Construct from root Region and set of subdivision ratios. @@ -111,6 +131,9 @@ class Subpaving return invalid; } + HostMirror + makeMirrorCopy() const; + /*! * @brief Get root region */ @@ -127,79 +150,53 @@ class Subpaving getDeviceMemorySize() const noexcept; /*! - * @brief Get tiles DualView - * @todo Rename this + * @brief Get the set of tiles */ - TilesDualView - getTilesView() + KOKKOS_INLINE_FUNCTION + const TilesView& + getTiles() const { return _tiles; } /*! - * @brief Synchronize Tile data onto specified memory space - */ - template - void - syncTiles(TContext context = onHost) - { - detail::syncUpdate(_tiles, context); - } - - /*! - * @brief Synchronize Zone data onto specified memory space - */ - template - void - syncZones(TContext context = onHost) - { - detail::syncUpdate(_zones, context); - } - - /*! - * @brief Synchronize Zone and Tile data onto specified memory space + * @brief Get random-access tiles view */ - template - void - syncAll(TContext context = onHost) + KOKKOS_INLINE_FUNCTION + const TilesRAView& + getRandomAccessTiles() const { - syncTiles(context); - syncZones(context); + return _tilesRA; } /*! - * @brief Get current number of tiles held in the specified memory space - * @note This does not synchronize + * @brief Get current number of tiles */ - template + KOKKOS_INLINE_FUNCTION IdType - getNumberOfTiles(TContext context = onHost) + getNumberOfTiles() const { - return static_cast(getTiles(context).extent(0)); + return static_cast(_tiles.size()); } /*! - * @brief Get the set of tiles in the specified memory space - * @note This does not synchronize + * @brief Get the set of zones */ - template KOKKOS_INLINE_FUNCTION - const TilesView& - getTiles(TContext context = onHost) const + const ZonesView& + getZones() const { - return detail::getContextualView(_tiles, context); + return _zones; } /*! - * @brief Get the set of zones in the specified memory space - * @note This does not synchronize + * @brief Get random-access zones view */ - template KOKKOS_INLINE_FUNCTION - const ZonesView& - getZones(TContext context = onHost) const + const ZonesRAView& + getRandomAccessZones() const { - return detail::getContextualView(_zones, context); + return _zonesRA; } /*! @@ -222,17 +219,31 @@ class Subpaving * @brief Perform a tree search (using the zones) for the given point, and * return the id of the containing tile (or invalid if not found) */ - template KOKKOS_INLINE_FUNCTION IdType - findTileId(const PointType& point, TContext context = onHost) const; + findTileId(const PointType& point) const; - /*! @cond */ +private: void - plot(); - /*! @endcond */ + setZones(const ZonesView& zones) + { + _zones = zones; + _zonesRA = _zones; + } + + void + setTiles(const TilesView& tiles) + { + _tiles = tiles; + _tilesRA = _tiles; + } + + void + setRefinementDepth(std::size_t depth) + { + _refinementDepth = depth; + } -private: /*! * @brief Check subdivision ratios for domain divisibility and copy final * form (one per level) into device view @@ -242,16 +253,35 @@ class Subpaving private: //! Zones represent the entire subdivision tree for the root region - ZonesDualView _zones; + ZonesView _zones; + ZonesRAView _zonesRA; //! Tiles represent the (selected) leaf nodes of the tree - TilesDualView _tiles; + TilesView _tiles; + TilesRAView _tilesRA; //! Region which fully encloses the domain of interest RegionType _rootRegion; //! Collection of SubdivisionInfo, one per expected refinement level - Kokkos::DualView*> _subdivisionInfos; + Kokkos::View*, MemorySpace> _subdivisionInfos; //! Level limit std::size_t _refinementDepth{}; }; + +namespace detail +{ +template +struct MemSpaceSubpavingHelper; + +template +struct MemSpaceSubpavingHelper> +{ + using Type = Subpaving; +}; +} // namespace detail + +template +using MemSpaceSubpaving = + typename detail::MemSpaceSubpavingHelper::Type; } // namespace plsm #include diff --git a/include/plsm/Subpaving.inl b/include/plsm/Subpaving.inl index 271b78f..cece801 100644 --- a/include/plsm/Subpaving.inl +++ b/include/plsm/Subpaving.inl @@ -12,26 +12,32 @@ namespace plsm { -template -Subpaving::Subpaving(const RegionType& region, +template +Subpaving::Subpaving( + const RegionType& region, const std::vector>& subdivisionRatios) : - _zones("zones", 1), _tiles("tiles", 1), _rootRegion(region) + _zones("zones", 1), + _zonesRA(_zones), + _tiles("tiles", 1), + _tilesRA(_tiles), + _rootRegion(region) { processSubdivisionRatios(subdivisionRatios); - _zones.h_view(0) = ZoneType{_rootRegion, 0}; - _zones.modify_host(); + auto zonesMirror = create_mirror_view(_zones); + zonesMirror[0] = ZoneType{_rootRegion, 0}; + deep_copy(_zones, zonesMirror); - _tiles.h_view(0) = TileType{_rootRegion, 0}; - _tiles.modify_host(); - - _zones.sync_device(); - _tiles.sync_device(); + auto tilesMirror = create_mirror_view(_tiles); + tilesMirror[0] = TileType{_rootRegion, 0}; + deep_copy(_tiles, tilesMirror); } -template +template void -Subpaving::processSubdivisionRatios( +Subpaving::processSubdivisionRatios( const std::vector>& subdivRatios) { auto subdivisionRatios = subdivRatios; @@ -44,23 +50,69 @@ Subpaving::processSubdivisionRatios( } return ret; }; - auto ratioProduct = std::accumulate(next(begin(subdivisionRatios)), - end(subdivisionRatios), subdivisionRatios.front(), elementWiseProduct); - auto getIntervalLength = [](const IntervalType& ival) { - return ival.length(); - }; + // Compute per-dimension ratio product for provided subdivision ratios + auto ratioProduct = + std::accumulate(begin(subdivisionRatios), end(subdivisionRatios), + SubdivisionRatio::filled(1), elementWiseProduct); + + // Get root region extents (what is being subdivided) std::array extents; std::transform(begin(_rootRegion), end(_rootRegion), begin(extents), - getIntervalLength); + [](auto&& ival) { return ival.length(); }); + + auto nonSelfFactors = [](auto x) { + using T = std::remove_reference_t; + std::vector result; + + // This will loop from 2 to sqrt(x) + for (T i = 2; i * i <= x; ++i) { + // Check if i divides x without leaving a remainder + if (x % i == 0) { + result.push_back(i); + // Include other factor if not root + if (x / i != i) { + result.push_back(x / i); + } + } + } + + std::sort(begin(result), end(result)); + return result; + }; + + auto getNextFactor = [nonSelfFactors](auto toSub, auto refFactor) { + using T = std::remove_reference_t; + if (toSub % refFactor == 0) { + return refFactor; + } + auto opt = nonSelfFactors(toSub); + if (opt.empty()) { + return static_cast(toSub); + } + T ret = 0; + for (auto i : opt) { + if (i > refFactor) { + continue; + } + ret = i; + } + if (ret == 0) { + ret = opt.front(); + } + return ret; + }; - // FIXME: infinite loop when ratios do not evenly divide extents + // Create additional ratio(s) until space is fully subdivided for (;;) { bool needAnotherLevel = false; - auto newRatio = SubdivisionRatio::filled(1); + auto nextRatio = SubdivisionRatio::filled(1); for (auto i : makeIntervalRange(Dim)) { if (ratioProduct[i] < extents[i]) { - newRatio[i] = subdivisionRatios.back()[i]; + // Determine next ratio to use to subdivide leftovers. Use last + // ratio if possible + nextRatio[i] = getNextFactor( + extents[i] / ratioProduct[i], subdivisionRatios.back()[i]); needAnotherLevel = true; } } @@ -69,8 +121,8 @@ Subpaving::processSubdivisionRatios( break; } - subdivisionRatios.push_back(newRatio); - ratioProduct = elementWiseProduct(ratioProduct, newRatio); + subdivisionRatios.push_back(nextRatio); + ratioProduct = elementWiseProduct(ratioProduct, nextRatio); } for (auto i : makeIntervalRange(Dim)) { @@ -83,34 +135,61 @@ Subpaving::processSubdivisionRatios( } } - _subdivisionInfos = Kokkos::DualView*>{ - "Subdivision Infos", subdivisionRatios.size()}; + _subdivisionInfos = + Kokkos::View*, MemorySpace>{ + "Subdivision Infos", subdivisionRatios.size()}; + auto subdivInfoMirror = create_mirror_view(_subdivisionInfos); std::copy(begin(subdivisionRatios), end(subdivisionRatios), - _subdivisionInfos.h_view.data()); - _subdivisionInfos.modify_host(); - _subdivisionInfos.sync_device(); + subdivInfoMirror.data()); + deep_copy(_subdivisionInfos, subdivInfoMirror); } -template +template +typename Subpaving::HostMirror +Subpaving::makeMirrorCopy() const +{ + HostMirror ret{}; + auto zones = create_mirror_view(_zones); + deep_copy(zones, _zones); + ret.setZones(zones); + + auto tiles = create_mirror_view(_tiles); + deep_copy(tiles, _tiles); + ret.setTiles(tiles); + + ret._rootRegion = _rootRegion; + + resize(ret._subdivisionInfos, _subdivisionInfos.size()); + deep_copy(ret._subdivisionInfos, _subdivisionInfos); + + ret._refinementDepth = _refinementDepth; + + return ret; +} + +template std::uint64_t -Subpaving::getDeviceMemorySize() const noexcept +Subpaving::getDeviceMemorySize() + const noexcept { std::uint64_t ret{}; - ret += _tiles.d_view.required_allocation_size(_tiles.d_view.extent(0)); - ret += _zones.d_view.required_allocation_size(_zones.d_view.extent(0)); + ret += _tiles.required_allocation_size(_tiles.size()); + ret += _zones.required_allocation_size(_zones.size()); ret += sizeof(_rootRegion); - ret += _subdivisionInfos.d_view.required_allocation_size( - _subdivisionInfos.d_view.extent(0)); + ret += _subdivisionInfos.required_allocation_size(_subdivisionInfos.size()); ret += sizeof(_refinementDepth); return ret; } -template +template template void -Subpaving::refine( +Subpaving::refine( TRefinementDetector&& detector) { using Refiner = detail::Refiner; @@ -118,16 +197,15 @@ Subpaving::refine( refiner(); } -template -template +template KOKKOS_INLINE_FUNCTION IdType -Subpaving::findTileId( - const PointType& point, TContext context) const +Subpaving::findTileId( + const PointType& point) const { - auto zones = getZones(context); IdType zoneId = 0; - auto zone = zones(zoneId); + auto zone = _zonesRA(zoneId); auto tileId = invalid; if (!zone.getRegion().contains(point)) { return tileId; @@ -139,7 +217,7 @@ Subpaving::findTileId( } auto newZoneId = zoneId; for (auto subZoneId : zone.getSubZoneRange()) { - if (zones(subZoneId).getRegion().contains(point)) { + if (_zonesRA(subZoneId).getRegion().contains(point)) { newZoneId = subZoneId; break; } @@ -148,37 +226,8 @@ Subpaving::findTileId( break; } zoneId = newZoneId; - zone = zones(zoneId); + zone = _zonesRA(zoneId); } return tileId; } - -//! @cond -template -void -Subpaving::plot() -{ - auto tiles = getTiles(); - std::ofstream ofs("gp.txt"); - for (auto i : makeIntervalRange(tiles.extent(0))) { - const auto& region = tiles(i).getRegion(); - ofs << "\n"; - ofs << region[0].begin() << " " << region[1].begin() << "\n"; - ofs << region[0].end() << " " << region[1].begin() << "\n"; - ofs << region[0].end() << " " << region[1].end() << "\n"; - ofs << region[0].begin() << " " << region[1].end() << "\n"; - ofs << region[0].begin() << " " << region[1].begin() << "\n"; - ofs << "\n"; - double q01 = 0.25 * region[0].begin() + 0.75 * region[0].end(); - double q03 = 0.75 * region[0].begin() + 0.25 * region[0].end(); - double q11 = 0.25 * region[1].begin() + 0.75 * region[1].end(); - double q13 = 0.75 * region[1].begin() + 0.25 * region[1].end(); - ofs << q01 << " " << q11 << "\n"; - ofs << q03 << " " << q13 << "\n"; - ofs << "\n"; - ofs << q01 << " " << q13 << "\n"; - ofs << q03 << " " << q11 << "\n"; - } -} -//! @endcond } // namespace plsm diff --git a/include/plsm/Utility.h b/include/plsm/Utility.h index 92f6483..4210202 100644 --- a/include/plsm/Utility.h +++ b/include/plsm/Utility.h @@ -4,7 +4,7 @@ #include #include -#include +#include #include @@ -20,11 +20,14 @@ inline constexpr T invalid = std::numeric_limits::max() - static_cast(1); template inline constexpr T wildcard = std::numeric_limits::max(); -template +template class Subpaving; -namespace detail -{ +using DefaultExecSpace = Kokkos::DefaultExecutionSpace; +using DefaultMemSpace = typename DefaultExecSpace::memory_space; +using DeviceMemSpace = DefaultMemSpace; +using HostMemSpace = Kokkos::HostSpace; + /*! * @brief Checks whether T is an instantiation of Subpaving */ @@ -38,13 +41,16 @@ struct IsSubpaving : std::false_type }; template -struct IsSubpaving<::plsm::Subpaving> : + typename TItemData, typename TMemSpace> +struct IsSubpaving< + ::plsm::Subpaving> : std::true_type { }; /*! @endcond */ +namespace detail +{ /*! * @brief Determine the type for the result of a subtraction */ diff --git a/include/plsm/Zone.h b/include/plsm/Zone.h index 0a1d168..7c6979c 100644 --- a/include/plsm/Zone.h +++ b/include/plsm/Zone.h @@ -41,6 +41,16 @@ struct Zone { } + /*! + * @brief Dimension of lattice + */ + static KOKKOS_INLINE_FUNCTION + constexpr DimType + dimension() noexcept + { + return RegionType::dimension(); + } + /*! * @brief Check if Zone owns a Tile (has a valid Tile index) */ diff --git a/include/plsm/detail/Refiner.h b/include/plsm/detail/Refiner.h index cd39fbd..d256449 100644 --- a/include/plsm/detail/Refiner.h +++ b/include/plsm/detail/Refiner.h @@ -4,9 +4,7 @@ #include #include -#include -#include #include #include #include @@ -16,6 +14,64 @@ namespace plsm { namespace detail { +struct ItemTotals +{ + IdType zones{0}; + IdType tiles{0}; + + KOKKOS_INLINE_FUNCTION + volatile ItemTotals& + operator+=(const volatile ItemTotals& other) volatile + { + zones += other.zones; + tiles += other.tiles; + return *this; + } +}; + +template +struct RefinerData +{ + static_assert(IsSubpaving{}); + + using SubpavingType = TSubpaving; + using ZoneType = typename SubpavingType::ZoneType; + using TileType = typename SubpavingType::TileType; + using ZonesView = typename SubpavingType::ZonesView; + using ZonesRAView = typename SubpavingType::ZonesRAView; + using TilesView = typename SubpavingType::TilesView; + + static constexpr DimType subpavingDim = SubpavingType::dimension(); + + using SubdivisionInfoType = SubdivisionInfo; + + using DetectorType = TDetector; + + ZonesView zones; + ZonesRAView zonesRA; + TilesView tiles; + + Kokkos::View subdivisionInfos; + + DetectorType detector; + + std::size_t targetDepth; + std::size_t currLevel{}; + + Kokkos::Array, subpavingDim> + enableRefine{}; + + Kokkos::View selectedSubZones{}; + + Kokkos::View newZoneCounts{}; + Kokkos::View subZoneStarts{}; + Kokkos::View newTileStarts{}; + + ItemTotals newItemTotals{}; + IdType numZones{zones.size()}; + IdType numTiles{tiles.size()}; +}; + /*! * @brief Refiner handles the refinement and selection of the Subpaving tiles */ @@ -23,7 +79,6 @@ template class Refiner { public: - static_assert(IsSubpaving::value, ""); using SubpavingType = TSubpaving; using ScalarType = typename SubpavingType::ScalarType; using ZoneType = typename SubpavingType::ZoneType; @@ -41,83 +96,26 @@ class Refiner void operator()(); - struct ItemTotals - { - IdType zones{0}; - IdType tiles{0}; - - KOKKOS_INLINE_FUNCTION - volatile ItemTotals& - operator+=(const volatile ItemTotals& other) volatile - { - zones += other.zones; - tiles += other.tiles; - return *this; - } - }; - - KOKKOS_INLINE_FUNCTION - IdType - countSelectSubZones(IdType index, const ZoneType& zone) const; - - KOKKOS_INLINE_FUNCTION void - countSelectNewItemsFromTile(IdType index, ItemTotals& runningTotals) const; - - void - countNewZonesAndTiles(); + countNewItems(); void findNewItemIndices(); - KOKKOS_INLINE_FUNCTION void - refineTile(IdType index) const; - - void - assignNewZonesAndTiles(); - - KOKKOS_INLINE_FUNCTION - RegionType - getSubZoneRegion(const ZoneType& zone, IdType subZoneLocalId, - const SubdivisionInfoType& subdivInfo) const; - - KOKKOS_INLINE_FUNCTION - SubdivisionRatioType - getSubdivisionRatio(std::size_t level, IdType tileIndex) const; + assignNewItems(); protected: - template + template friend class ::plsm::Subpaving; Refiner(SubpavingType& subpaving, const DetectorType& detector); protected: SubpavingType& _subpaving; - using ZonesView = typename SubpavingType::template ZonesView; - ZonesView _zones; - using TilesView = typename SubpavingType::template TilesView; - TilesView _tiles; - - Kokkos::DualView _subdivisionInfos; - - std::size_t _currLevel; - std::size_t _targetDepth; - - DetectorType _detector; - - using DefaultExecSpace = - typename Kokkos::View::traits::execution_space; - Kokkos::Array, subpavingDim> _enableRefine; - - Kokkos::View _selectedSubZones; + typename Kokkos::View::HostMirror _subdivInfoMirror; - ItemTotals _newItemTotals{}; - IdType _numTiles; - IdType _numZones; - Kokkos::View _newZoneCounts; - Kokkos::View _subZoneStarts; - Kokkos::View _newTileStarts; + RefinerData _data; }; } // namespace detail } // namespace plsm diff --git a/include/plsm/detail/Refiner.inl b/include/plsm/detail/Refiner.inl index 0e2710b..4245310 100644 --- a/include/plsm/detail/Refiner.inl +++ b/include/plsm/detail/Refiner.inl @@ -8,98 +8,62 @@ namespace plsm { namespace detail { -template -class RefinerFunctor : public TRefiner -{ -public: - RefinerFunctor(const TRefiner& refiner) : TRefiner(refiner) - { - } -}; - -template -class CountNewItemsFromTile : public RefinerFunctor -{ -public: - using ItemTotals = typename TRefiner::ItemTotals; - - using RefinerFunctor::RefinerFunctor; +using AllocNoInit = Kokkos::ViewAllocateWithoutInitializing; - KOKKOS_INLINE_FUNCTION - void - operator()(IdType index, ItemTotals& running) const - { - this->countSelectNewItemsFromTile(index, running); - } -}; - -template -class RefineTile : public RefinerFunctor +template +KOKKOS_INLINE_FUNCTION +SubdivisionRatio +getSubdivisionRatio(const TData& data, std::size_t level, IdType tileIndex) { -public: - using RefinerFunctor::RefinerFunctor; - - KOKKOS_INLINE_FUNCTION - void - operator()(IdType index) const - { - this->refineTile(index); + auto ret = data.subdivisionInfos[level].getRatio(); + for (DimType i = 0; i < data.subpavingDim; ++i) { + if (!data.enableRefine[i].test(static_cast(tileIndex))) { + ret[i] = 1; + } } -}; - -template -Refiner::Refiner( - SubpavingType& subpaving, const DetectorType& detector) : - _subpaving(subpaving), - _zones(subpaving._zones.d_view), - _tiles(subpaving._tiles.d_view), - _subdivisionInfos(subpaving._subdivisionInfos), - _detector(detector), - _numTiles(static_cast(_tiles.extent(0))), - _numZones(static_cast(_zones.extent(0))) -{ + return ret; } -template -void -Refiner::operator()() +template +KOKKOS_INLINE_FUNCTION +typename TZone::RegionType +getSubZoneRegion(const TZone& zone, IdType subZoneLocalId, + const SubdivisionInfo& subdivInfo) { - _targetDepth = (_detector.depth() == _detector.fullDepth) ? - this->_subdivisionInfos.h_view.extent(0) : - _detector.depth(); - - for (_currLevel = 0; _currLevel < _targetDepth; ++_currLevel) { - countNewZonesAndTiles(); + using RegionType = typename TZone::RegionType; + using ScalarType = typename RegionType::ScalarType; + using IntervalType = typename RegionType::IntervalType; - if (_newItemTotals.zones == 0) { - break; - } + constexpr auto subpavingDim = TZone::dimension(); - findNewItemIndices(); + MultiIndex mId = subdivInfo.getMultiIndex(subZoneLocalId); - assignNewZonesAndTiles(); + const auto& zoneRegion = zone.getRegion(); + RegionType ret; + for (auto i : makeIntervalRange(subpavingDim)) { + const auto& ival = zoneRegion[i]; + auto delta = ival.length() / subdivInfo.getRatio()[i]; + ret[i] = + IntervalType{ival.begin() + static_cast(mId[i] * delta), + ival.begin() + static_cast((mId[i] + 1) * delta)}; } - - _subpaving._zones.d_view = _zones; - _subpaving._zones.modify_device(); - _subpaving._tiles.d_view = _tiles; - _subpaving._tiles.modify_device(); - _subpaving._refinementDepth = _currLevel; + return ret; } -template +template KOKKOS_INLINE_FUNCTION IdType -Refiner::countSelectSubZones( - IdType index, const ZoneType& zone) const +countSelectSubZones( + const TData& data, IdType index, const typename TData::ZoneType& zone) { - SubdivisionInfoType info(getSubdivisionRatio(zone.getLevel(), index)); + auto info = SubdivisionInfo{ + getSubdivisionRatio(data, zone.getLevel(), index)}; auto numSubRegions = info.getRatio().getProduct(); - auto selected = Kokkos::subview(_selectedSubZones, index, Kokkos::ALL); + auto selected = Kokkos::subview(data.selectedSubZones, index, Kokkos::ALL); IdType count = 0; for (auto i : makeIntervalRange(numSubRegions)) { auto subRegion = getSubZoneRegion(zone, i, info); - if (_detector(DetectorType::selectTag, subRegion)) { + if (data.detector(data.detector.selectTag, subRegion)) { selected(count) = i; ++count; } @@ -107,85 +71,171 @@ Refiner::countSelectSubZones( return count; } -template +template KOKKOS_INLINE_FUNCTION void -Refiner::countSelectNewItemsFromTile( - IdType index, ItemTotals& runningTotals) const +countNewItemsFromTile( + const TData& data, IdType index, ItemTotals& runningTotals) { - using BoolVec = typename DetectorType::template BoolVec; - - const auto& tile = _tiles(index); + using RegionType = typename TData::ZoneType::RegionType; + using BoolVec = refine::BoolVec; + const auto& tile = data.tiles(index); auto zoneId = tile.getOwningZoneIndex(); IdType count = 0; - auto& zone = _zones(zoneId); + auto& zone = data.zonesRA(zoneId); auto level = zone.getLevel(); - if (level < _targetDepth) { + if (level < data.targetDepth) { BoolVec enable{}; - if (_detector(DetectorType::refineTag, tile.getRegion(), enable)) { - for (DimType i = 0; i < subpavingDim; ++i) { + if (data.detector(data.detector.refineTag, tile.getRegion(), enable)) { + for (DimType i = 0; i < data.subpavingDim; ++i) { if (enable[i]) { - _enableRefine[i].set(static_cast(index)); + data.enableRefine[i].set(static_cast(index)); } else { - _enableRefine[i].reset(static_cast(index)); + data.enableRefine[i].reset(static_cast(index)); } } - count = countSelectSubZones(index, zone); + count = countSelectSubZones(data, index, zone); } } - _newZoneCounts(index) = count; + data.newZoneCounts(index) = count; if (count > 0) { runningTotals.zones += count; runningTotals.tiles += count - 1; } } +template +KOKKOS_INLINE_FUNCTION +void +refineTile(const TData& data, IdType index) +{ + using ZoneType = typename TData::ZoneType; + using TileType = typename TData::TileType; + + auto newZones = data.newZoneCounts(index); + if (newZones == 0) { + return; + } + + auto& tile = data.tiles(index); + auto ownerZoneId = tile.getOwningZoneIndex(); + auto& ownerZone = data.zones(ownerZoneId); + auto level = ownerZone.getLevel(); + auto newLevel = level + 1; + auto info = SubdivisionInfo{ + getSubdivisionRatio(data, level, index)}; + + // Create first new zone, replace current tile and associate + auto subZoneBeginId = data.numZones + data.subZoneStarts(index); + data.zones(subZoneBeginId) = ZoneType{ + getSubZoneRegion(ownerZone, data.selectedSubZones(index, 0), info), + newLevel, ownerZoneId}; + data.zones(subZoneBeginId).setTileIndex(index); + tile = TileType{data.zonesRA(subZoneBeginId).getRegion(), subZoneBeginId}; + + // Create and associate remaining zones and tiles + auto tileBeginId = data.numTiles + data.newTileStarts(index); + for (IdType i = 1; i < newZones; ++i) { + auto zoneId = subZoneBeginId + i; + auto tileId = tileBeginId + i - 1; + data.zones(zoneId) = ZoneType{ + getSubZoneRegion(ownerZone, data.selectedSubZones(index, i), info), + newLevel, ownerZoneId}; + data.zones(zoneId).setTileIndex(tileId); + data.tiles(tileId) = TileType{data.zonesRA(zoneId).getRegion(), zoneId}; + } + + ownerZone.removeTile(); + ownerZone.setSubZoneIndices({subZoneBeginId, subZoneBeginId + newZones}); +} + +template +Refiner::Refiner( + SubpavingType& subpaving, const DetectorType& detector) : + _subpaving(subpaving), + _subdivInfoMirror(create_mirror_view(subpaving._subdivisionInfos)), + _data{subpaving._zones, subpaving._zonesRA, subpaving._tiles, + subpaving._subdivisionInfos, detector, + (detector.depth() == detector.fullDepth) ? + subpaving._subdivisionInfos.size() : + detector.depth()} +{ + deep_copy(_subdivInfoMirror, _data.subdivisionInfos); +} + +template +void +Refiner::operator()() +{ + for (_data.currLevel = 0; _data.currLevel < _data.targetDepth; + ++_data.currLevel) { + countNewItems(); + + if (_data.newItemTotals.zones == 0) { + break; + } + + findNewItemIndices(); + + assignNewItems(); + } + + _subpaving.setZones(_data.zones); + _subpaving.setTiles(_data.tiles); + _subpaving.setRefinementDepth(_data.currLevel); +} + template void -Refiner::countNewZonesAndTiles() +Refiner::countNewItems() { - _newZoneCounts = Kokkos::View( - Kokkos::ViewAllocateWithoutInitializing{"New Zone Counts"}, _numTiles); + auto numTiles = _data.tiles.size(); + _data.newZoneCounts = + Kokkos::View(AllocNoInit{"New Zone Counts"}, numTiles); auto numSubZones = - _subdivisionInfos.h_view(_currLevel).getRatio().getProduct(); - _selectedSubZones = Kokkos::View( - Kokkos::ViewAllocateWithoutInitializing{"Selected Sub-Zones"}, - _numTiles, numSubZones); - std::for_each(begin(_enableRefine), end(_enableRefine), - [numTiles = _numTiles](auto&& bitset) { - bitset = Kokkos::Bitset( - static_cast(numTiles)); + _subdivInfoMirror[_data.currLevel].getRatio().getProduct(); + _data.selectedSubZones = Kokkos::View( + AllocNoInit{"Selected Sub-Zones"}, numTiles, numSubZones); + std::for_each(begin(_data.enableRefine), end(_data.enableRefine), + [numTiles](auto&& bitset) { + using Bitset = std::remove_reference_t; + bitset = Bitset(static_cast(numTiles)); }); ItemTotals counts{}; - Kokkos::parallel_reduce("CountNewItemsFromTile", - _numTiles, CountNewItemsFromTile{*this}, counts); + auto data = _data; + Kokkos::parallel_reduce( + "CountNewItemsFromTile", numTiles, + KOKKOS_LAMBDA(IdType id, ItemTotals & running) { + countNewItemsFromTile(data, id, running); + }, + counts); Kokkos::fence(); - _newItemTotals = counts; + _data.newItemTotals = counts; } template void Refiner::findNewItemIndices() { - auto subZoneStarts = Kokkos::View( - Kokkos::ViewAllocateWithoutInitializing{"SubZone Start Ids"}, - _numTiles); - auto newTileStarts = Kokkos::View( - Kokkos::ViewAllocateWithoutInitializing{"Tile Start Ids"}, _numTiles); - auto newZoneCounts = _newZoneCounts; + auto numTiles = _data.tiles.size(); + auto subZoneStarts = + Kokkos::View(AllocNoInit{"SubZone Start Ids"}, numTiles); + auto newTileStarts = + Kokkos::View(AllocNoInit{"Tile Start Ids"}, numTiles); + auto newZoneCounts = _data.newZoneCounts; // Initialize starts - Kokkos::parallel_for("InitializeNewItemStarts", - _numTiles, KOKKOS_LAMBDA(IdType i) { + Kokkos::parallel_for( + "InitializeNewItemStarts", numTiles, KOKKOS_LAMBDA(IdType i) { auto newZoneCount = newZoneCounts(i); subZoneStarts(i) = newZoneCount; newTileStarts(i) = (newZoneCount == 0) ? 0 : newZoneCount - 1; }); ItemTotals totals{}; - Kokkos::parallel_scan("ScanNewItemStarts", - _numTiles, + Kokkos::parallel_scan( + "ScanNewItemStarts", numTiles, KOKKOS_LAMBDA(IdType i, ItemTotals & update, const bool finalPass) { const auto tmpZones = subZoneStarts(i); const auto tmpTiles = newTileStarts(i); @@ -199,100 +249,26 @@ Refiner::findNewItemIndices() totals); Kokkos::fence(); - _subZoneStarts = subZoneStarts; - _newTileStarts = newTileStarts; -} - -template -KOKKOS_INLINE_FUNCTION -typename Refiner::RegionType -Refiner::getSubZoneRegion(const ZoneType& zone, - IdType subZoneLocalId, const SubdivisionInfoType& subdivInfo) const -{ - using IntervalType = typename RegionType::IntervalType; - - MultiIndex mId = subdivInfo.getMultiIndex(subZoneLocalId); - - const auto& zoneRegion = zone.getRegion(); - RegionType ret; - for (auto i : makeIntervalRange(subpavingDim)) { - const auto& ival = zoneRegion[i]; - auto delta = ival.length() / subdivInfo.getRatio()[i]; - ret[i] = - IntervalType{ival.begin() + static_cast(mId[i] * delta), - ival.begin() + static_cast((mId[i] + 1) * delta)}; - } - return ret; + _data.subZoneStarts = subZoneStarts; + _data.newTileStarts = newTileStarts; } template -KOKKOS_INLINE_FUNCTION -typename Refiner::SubdivisionRatioType -Refiner::getSubdivisionRatio( - std::size_t level, IdType tileIndex) const -{ - SubdivisionRatioType ret = _subdivisionInfos.d_view[level].getRatio(); - for (DimType i = 0; i < subpavingDim; ++i) { - if (!_enableRefine[i].test(static_cast(tileIndex))) { - ret[i] = 1; - } - } - return ret; -} - -template -KOKKOS_INLINE_FUNCTION void -Refiner::refineTile(IdType index) const +Refiner::assignNewItems() { - auto newZones = _newZoneCounts(index); - if (newZones == 0) { - return; - } - - auto& tile = _tiles(index); - auto ownerZoneId = tile.getOwningZoneIndex(); - auto& ownerZone = _zones(ownerZoneId); - auto level = ownerZone.getLevel(); - auto newLevel = level + 1; - auto info = SubdivisionInfoType(getSubdivisionRatio(level, index)); - - // Create first new zone, replace current tile and associate - auto subZoneBeginId = _numZones + _subZoneStarts(index); - _zones(subZoneBeginId) = - ZoneType{getSubZoneRegion(ownerZone, _selectedSubZones(index, 0), info), - newLevel, ownerZoneId}; - _zones(subZoneBeginId).setTileIndex(index); - tile = TileType{_zones(subZoneBeginId).getRegion(), subZoneBeginId}; - - // Create and associate remaining zones and tiles - auto tileBeginId = _numTiles + _newTileStarts(index); - for (IdType i = 1; i < newZones; ++i) { - auto zoneId = subZoneBeginId + i; - auto tileId = tileBeginId + i - 1; - _zones(zoneId) = ZoneType{ - getSubZoneRegion(ownerZone, _selectedSubZones(index, i), info), - newLevel, ownerZoneId}; - _zones(zoneId).setTileIndex(tileId); - _tiles(tileId) = TileType{_zones(zoneId).getRegion(), zoneId}; - } - - ownerZone.removeTile(); - ownerZone.setSubZoneIndices({subZoneBeginId, subZoneBeginId + newZones}); -} - -template -void -Refiner::assignNewZonesAndTiles() -{ - Kokkos::resize(_zones, _zones.extent(0) + _newItemTotals.zones); - Kokkos::resize(_tiles, _tiles.extent(0) + _newItemTotals.tiles); - - Kokkos::parallel_for("RefineTile", _numTiles, RefineTile{*this}); + Kokkos::resize(_data.zones, _data.numZones + _data.newItemTotals.zones); + _data.zonesRA = _data.zones; + Kokkos::resize(_data.tiles, _data.numTiles + _data.newItemTotals.tiles); + + auto data = _data; + Kokkos::parallel_for( + "RefineTile", data.numTiles, + KOKKOS_LAMBDA(IdType id) { refineTile(data, id); }); Kokkos::fence(); - _numTiles = static_cast(_tiles.extent(0)); - _numZones = static_cast(_zones.extent(0)); + _data.numZones = _data.zones.size(); + _data.numTiles = _data.tiles.size(); } } // namespace detail } // namespace plsm diff --git a/include/plsm/refine/Detector.h b/include/plsm/refine/Detector.h index aabddda..a97ef05 100644 --- a/include/plsm/refine/Detector.h +++ b/include/plsm/refine/Detector.h @@ -8,6 +8,12 @@ namespace plsm { namespace refine { +/*! + * @brief Set of bool flags, one for each axis + */ +template +using BoolVec = Kokkos::Array; + /*! * Tag for straightforward refinement decision */ @@ -169,12 +175,6 @@ class Detector //! Tag specifying selection method using SelectTag = typename Traits::SelectTag; - /*! - * @brief Set of bool flags, one for each axis - */ - template - using BoolVec = Kokkos::Array; - //! Tag instance for refinement static constexpr RefineTag refineTag{}; //! Tag instance for selection diff --git a/include/plsm/refine/MultiDetector.h b/include/plsm/refine/MultiDetector.h index e85d6c7..7fc8d5d 100644 --- a/include/plsm/refine/MultiDetector.h +++ b/include/plsm/refine/MultiDetector.h @@ -60,9 +60,6 @@ class MultiDetectorImpl : //! Alias for next link in the chain using Tail = MultiDetectorImpl; - template - using BoolVec = typename Head::template BoolVec; - MultiDetectorImpl(const MultiDetectorImpl&) = default; MultiDetectorImpl(MultiDetectorImpl&&) = default; @@ -128,9 +125,6 @@ class MultiDetector : public Detector> //! Alias for parent class type using Superclass = Detector>; - template - using BoolVec = typename Superclass::template BoolVec; - MultiDetector(const MultiDetector&) = default; MultiDetector(MultiDetector&&) = default; diff --git a/test/benchmark_Subpaving.cpp b/test/benchmark_Subpaving.cpp index 19d0961..acf0cd9 100644 --- a/test/benchmark_Subpaving.cpp +++ b/test/benchmark_Subpaving.cpp @@ -26,7 +26,6 @@ TEST_CASE("Subpaving 2D", "[Subpaving]") using Tags = TagPair; s.refine(BallDetector{{0, 0}, 500}); }; - std::cout << "depth: " << s.getRefinementDepth() << '\n'; test::renderSubpaving(s); } @@ -55,7 +54,6 @@ TEST_CASE("Subpaving 2D", "[Subpaving]") PD{{{{0, 0}}, {{256, 128}}, {{384, 256}}, {{512, 512}}}}, RD{{Ival{0, 56}, Ival{0, 56}}})); }; - std::cout << "z-aligned\n"; test::renderSubpaving(spv); } } @@ -78,6 +76,21 @@ TEST_CASE("Subpaving 3D", "[Subpaving]") s.refine(BallDetector{{256, 256, 256}, 128}); }; test::renderSubpaving(s); + std::size_t errors = 0; + BENCHMARK("search: ball") + { + auto tiles = s.getTiles(); + Kokkos::parallel_reduce( + tiles.size(), + KOKKOS_LAMBDA(std::size_t i, std::size_t & running) { + auto id = s.findTileId(tiles[i].getRegion().getOrigin()); + if (id != i) { + ++running; + } + }, + errors); + }; + REQUIRE(errors == 0); } SECTION("z-aligned") @@ -112,6 +125,21 @@ TEST_CASE("Subpaving 3D", "[Subpaving]") RegionDetector{{ival, ival, ival}})); }; test::renderSubpaving(s); + std::size_t errors = 0; + BENCHMARK("search: z-aligned polyline plus box") + { + auto tiles = s.getTiles(); + Kokkos::parallel_reduce( + tiles.size(), + KOKKOS_LAMBDA(std::size_t i, std::size_t & running) { + auto id = s.findTileId(tiles[i].getRegion().getOrigin()); + if (id != i) { + ++running; + } + }, + errors); + }; + REQUIRE(errors == 0); } SECTION("x-aligned") @@ -183,17 +211,32 @@ TEST_CASE("Subpaving with XRN Defaults", "[Subpaving][XRN]") { using RegionType = typename Subpaving::RegionType; using Ival = Interval; - RegionType r{{Ival{0, 5120}, Ival{0, 4096}, Ival{0, 16}}}; - Subpaving s(r, {{{10, 8, 2}}, {{8, 8, 2}}}); + RegionType r{{Ival{0, 10240}, Ival{0, 8192}, Ival{0, 16}}}; + Subpaving s(r, {{{10, 8, 2}}, {{2, 2, 2}}}); std::vector> rspecPoints; constexpr auto wild = wildcard; // rspecPoints.push_back({{wild, wild, 3}}); rspecPoints.push_back({{0, 0, wild}}); rspecPoints.push_back({{1000, 500, wild}}); - rspecPoints.push_back({{5120, 4096, wild}}); - BENCHMARK("refine: XRN Default") + rspecPoints.push_back({{10240, 8192, wild}}); + BENCHMARK("refine: XRN") { s.refine(refine::PolylineDetector{rspecPoints}); }; test::renderSubpaving(s); + std::size_t errors = 0; + BENCHMARK("search: XRN") + { + auto tiles = s.getTiles(); + Kokkos::parallel_reduce( + tiles.size(), + KOKKOS_LAMBDA(std::size_t i, std::size_t & running) { + auto id = s.findTileId(tiles[i].getRegion().getOrigin()); + if (id != i) { + ++running; + } + }, + errors); + }; + REQUIRE(errors == 0); } diff --git a/test/include/plsm/PrintSubpaving.h b/test/include/plsm/PrintSubpaving.h index 546db16..6d58be6 100644 --- a/test/include/plsm/PrintSubpaving.h +++ b/test/include/plsm/PrintSubpaving.h @@ -8,11 +8,11 @@ namespace plsm { namespace test { -template +template void -printSubpaving(Subpaving& subpaving, std::ostream& os) +printSubpaving(TSubpaving& subpaving, std::ostream& os) { - subpaving.syncAll(onHost); + static_assert(plsm::IsSubpaving{}); auto tiles = subpaving.getTiles(); auto zones = subpaving.getZones(); os << "Zones: " << zones.size() << '\n'; diff --git a/test/include/plsm/RenderSubpaving.h b/test/include/plsm/RenderSubpaving.h index b98401c..5ce7592 100644 --- a/test/include/plsm/RenderSubpaving.h +++ b/test/include/plsm/RenderSubpaving.h @@ -17,11 +17,12 @@ namespace plsm { namespace test { -template -void -renderSubpaving(Subpaving& subpaving) +template +inline void +renderSubpaving(Subpaving& sp) { - subpaving.syncTiles(onHost); + auto subpaving = sp.makeMirrorCopy(); auto tiles = subpaving.getTiles(); auto numTiles = tiles.extent(0); @@ -49,6 +50,47 @@ renderSubpaving(Subpaving& subpaving) region[0].end(), region[1].end(), region[2].end()); points->InsertNextPoint( region[0].begin(), region[1].end(), region[2].end()); + + cells->InsertNextCell(2); + cells->InsertCellPoint(pId); + cells->InsertCellPoint(pId + 1); + cells->InsertNextCell(2); + cells->InsertCellPoint(pId + 1); + cells->InsertCellPoint(pId + 2); + cells->InsertNextCell(2); + cells->InsertCellPoint(pId + 2); + cells->InsertCellPoint(pId + 3); + cells->InsertNextCell(2); + cells->InsertCellPoint(pId + 3); + cells->InsertCellPoint(pId); + + cells->InsertNextCell(2); + cells->InsertCellPoint(pId + 4); + cells->InsertCellPoint(pId + 5); + cells->InsertNextCell(2); + cells->InsertCellPoint(pId + 5); + cells->InsertCellPoint(pId + 6); + cells->InsertNextCell(2); + cells->InsertCellPoint(pId + 6); + cells->InsertCellPoint(pId + 7); + cells->InsertNextCell(2); + cells->InsertCellPoint(pId + 7); + cells->InsertCellPoint(pId + 4); + + cells->InsertNextCell(2); + cells->InsertCellPoint(pId + 0); + cells->InsertCellPoint(pId + 4); + cells->InsertNextCell(2); + cells->InsertCellPoint(pId + 1); + cells->InsertCellPoint(pId + 5); + cells->InsertNextCell(2); + cells->InsertCellPoint(pId + 2); + cells->InsertCellPoint(pId + 6); + cells->InsertNextCell(2); + cells->InsertCellPoint(pId + 3); + cells->InsertCellPoint(pId + 7); + + pId += 8; } else if (Dim == 2) { points->InsertNextPoint(region[0].begin(), region[1].begin(), 0.0); @@ -56,55 +98,24 @@ renderSubpaving(Subpaving& subpaving) points->InsertNextPoint(region[0].end(), region[1].end(), 0.0); points->InsertNextPoint(region[0].begin(), region[1].end(), 0.0); - points->InsertNextPoint(region[0].begin(), region[1].begin(), 0.0); - points->InsertNextPoint(region[0].end(), region[1].begin(), 0.0); - points->InsertNextPoint(region[0].end(), region[1].end(), 0.0); - points->InsertNextPoint(region[0].begin(), region[1].end(), 0.0); + cells->InsertNextCell(2); + cells->InsertCellPoint(pId); + cells->InsertCellPoint(pId + 1); + cells->InsertNextCell(2); + cells->InsertCellPoint(pId + 1); + cells->InsertCellPoint(pId + 2); + cells->InsertNextCell(2); + cells->InsertCellPoint(pId + 2); + cells->InsertCellPoint(pId + 3); + cells->InsertNextCell(2); + cells->InsertCellPoint(pId + 3); + cells->InsertCellPoint(pId); + + pId += 4; } else { throw std::runtime_error("Unsupported dimensionality"); } - - cells->InsertNextCell(2); - cells->InsertCellPoint(pId); - cells->InsertCellPoint(pId + 1); - cells->InsertNextCell(2); - cells->InsertCellPoint(pId + 1); - cells->InsertCellPoint(pId + 2); - cells->InsertNextCell(2); - cells->InsertCellPoint(pId + 2); - cells->InsertCellPoint(pId + 3); - cells->InsertNextCell(2); - cells->InsertCellPoint(pId + 3); - cells->InsertCellPoint(pId); - - cells->InsertNextCell(2); - cells->InsertCellPoint(pId + 4); - cells->InsertCellPoint(pId + 5); - cells->InsertNextCell(2); - cells->InsertCellPoint(pId + 5); - cells->InsertCellPoint(pId + 6); - cells->InsertNextCell(2); - cells->InsertCellPoint(pId + 6); - cells->InsertCellPoint(pId + 7); - cells->InsertNextCell(2); - cells->InsertCellPoint(pId + 7); - cells->InsertCellPoint(pId + 4); - - cells->InsertNextCell(2); - cells->InsertCellPoint(pId + 0); - cells->InsertCellPoint(pId + 4); - cells->InsertNextCell(2); - cells->InsertCellPoint(pId + 1); - cells->InsertCellPoint(pId + 5); - cells->InsertNextCell(2); - cells->InsertCellPoint(pId + 2); - cells->InsertCellPoint(pId + 6); - cells->InsertNextCell(2); - cells->InsertCellPoint(pId + 3); - cells->InsertCellPoint(pId + 7); - - pId += 8; } auto pd = vtkSmartPointer::New(); @@ -123,9 +134,9 @@ renderSubpaving(Subpaving& subpaving) static auto interactor = [=]() { auto window = vtkSmartPointer::New(); window->AddRenderer(renderer); - auto interactor = vtkSmartPointer::New(); - interactor->SetRenderWindow(window); - return interactor; + auto ia = vtkSmartPointer::New(); + ia->SetRenderWindow(window); + return ia; }(); interactor->GetRenderWindow()->Render(); interactor->Start(); @@ -139,14 +150,49 @@ namespace plsm { namespace test { -template -void -renderSubpaving(Subpaving& subpaving) +template +inline void +renderSubpaving(Subpaving& subpaving) { - std::cout << "Number of Tiles: " << subpaving.getNumberOfTiles(onDevice) + std::cout << "\nNumber of Tiles: " << subpaving.getNumberOfTiles() << std::endl; } } // namespace test } // namespace plsm #endif + +namespace plsm +{ +namespace test +{ +template +inline void +plot(Subpaving& subpaving) +{ + auto tiles = subpaving.getTiles(); + std::ofstream ofs("gp.txt"); + for (auto i : makeIntervalRange(tiles.extent(0))) { + const auto& region = tiles(i).getRegion(); + ofs << "\n"; + ofs << region[0].begin() << " " << region[1].begin() << "\n"; + ofs << region[0].end() << " " << region[1].begin() << "\n"; + ofs << region[0].end() << " " << region[1].end() << "\n"; + ofs << region[0].begin() << " " << region[1].end() << "\n"; + ofs << region[0].begin() << " " << region[1].begin() << "\n"; + ofs << "\n"; + double q01 = 0.25 * region[0].begin() + 0.75 * region[0].end(); + double q03 = 0.75 * region[0].begin() + 0.25 * region[0].end(); + double q11 = 0.25 * region[1].begin() + 0.75 * region[1].end(); + double q13 = 0.75 * region[1].begin() + 0.25 * region[1].end(); + ofs << q01 << " " << q11 << "\n"; + ofs << q03 << " " << q13 << "\n"; + ofs << "\n"; + ofs << q01 << " " << q13 << "\n"; + ofs << q03 << " " << q11 << "\n"; + } +} +} // namespace test +} // namespace plsm diff --git a/test/unittest_Detectors.cpp b/test/unittest_Detectors.cpp index 8d8191b..a214160 100644 --- a/test/unittest_Detectors.cpp +++ b/test/unittest_Detectors.cpp @@ -407,7 +407,7 @@ TEMPLATE_LIST_TEST_CASE( Kokkos::parallel_reduce( 1, KOKKOS_LAMBDA(std::size_t, int& fl) { - typename PD::template BoolVec> res; + BoolVec> res; K_REQUIRE(md.refine(r1, res), fl); K_REQUIRE(md.refine(r2, res), fl); K_REQUIRE(!md.refine(r3, res), fl); diff --git a/test/unittest_Subpaving.cpp b/test/unittest_Subpaving.cpp index 040c632..20e3798 100644 --- a/test/unittest_Subpaving.cpp +++ b/test/unittest_Subpaving.cpp @@ -10,6 +10,54 @@ #include using namespace plsm; +namespace plsm::test +{ +template +struct SubpavingTester +{ + using Ratio = SubdivisionRatio; + + void + checkRatios(const std::vector& ratios) + { + auto infos = subpaving.makeMirrorCopy()._subdivisionInfos; + REQUIRE(infos.size() == ratios.size()); + } + + TSubpaving subpaving; +}; + +template +SubpavingTester(const TSubpaving&) -> SubpavingTester; +} // namespace plsm::test + +TEMPLATE_LIST_TEST_CASE( + "Process Subdivision Ratios", "[Subpaving][template]", test::IntTypes) +{ + auto subpaving = Subpaving({{{0, 100}, {0, 100}}}, {{{5, 5}}}); + test::SubpavingTester{subpaving}.checkRatios( + {{{5, 5}, {5, 5}, {2, 2}, {2, 2}}}); + + subpaving = Subpaving({{{0, 250}, {0, 50}}}, {{{5, 5}}}); + test::SubpavingTester{subpaving}.checkRatios( + {{{5, 5}, {5, 5}, {5, 2}, {2, 2}}}); + + subpaving = Subpaving({{{10, 20}, {25, 35}}}, {{{5, 5}}}); + test::SubpavingTester{subpaving}.checkRatios({{{5, 5}, {2, 2}}}); + + subpaving = Subpaving({{{0, 300}, {0, 275}}}, {{{5, 5}}}); + test::SubpavingTester{subpaving}.checkRatios( + {{{5, 5}, {5, 5}, {4, 11}, {3, 1}}}); + + subpaving = Subpaving({{{0, 2116}, {0, 1155}}}, {{{23, 11}}}); + test::SubpavingTester{subpaving}.checkRatios( + {{{23, 11}, {23, 7}, {2, 5}, {2, 3}}}); + + subpaving = Subpaving({{{0, 2116}, {0, 1155}}}, {{{2, 3}}}); + test::SubpavingTester{subpaving}.checkRatios( + {{{2, 3}, {2, 5}, {23, 7}, {23, 11}}}); +} + TEMPLATE_LIST_TEST_CASE( "Subpaving Basic", "[Subpaving][template]", test::IntTypes) { @@ -24,22 +72,20 @@ TEMPLATE_LIST_TEST_CASE( using RegionDetector = refine::RegionDetector>; sp.refine(RegionDetector{sp.getLatticeRegion()}); - REQUIRE(sp.getTiles(onDevice).extent(0) == 64); - sp.syncTiles(onHost); REQUIRE(sp.getTiles().extent(0) == 64); - REQUIRE(sp.findTileId({3, 3, 3}) == invalid); - sp.syncAll(onHost); - REQUIRE(sp.findTileId({3, 3, 3}) == 63); + + auto sph = sp.makeMirrorCopy(); + REQUIRE(sph.findTileId({3, 3, 3}) == 63); using Range3D = Kokkos::MDRangePolicy>; std::size_t errors = 0; - auto tiles = sp.getTiles(onDevice); + auto tiles = sp.getTiles(); Kokkos::parallel_reduce( Range3D({0, 0, 0}, {4, 4, 4}), KOKKOS_LAMBDA( TestType i, TestType j, TestType k, std::size_t & running) { - auto tileId = sp.findTileId(PointType{i, j, k}, onDevice); + auto tileId = sp.findTileId(PointType{i, j, k}); if (tileId == invalid) { ++running; } @@ -53,7 +99,7 @@ TEMPLATE_LIST_TEST_CASE( KOKKOS_LAMBDA( TestType i, TestType j, TestType k, std::size_t & running) { PointType p({i, j, k}); - auto tileId = sp.findTileId(p, onDevice); + auto tileId = sp.findTileId(p); if (tiles(tileId).getRegion().getOrigin() != p) { ++running; }